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ABSTRACT 


We investigated the influence of transport phenomena on the 
morphology of crystalline materials. Two problems were studied: 
one dealt with the effects of convection on the crystallization 
of pure materials, the other with the crystal 1 i z at i on of proteins 
from solution. In the first study we were interested in how 
convection alters the stability of the growth process and the 
relation between undercooling and the growth speed. In the 
second, we sought to find out why protein crystals grow as slowly 
as they do and how crystal morphology depends on the growth rate 
and crystal size. In both studies the research focussed on 
developing fundamental data that is a prerequisite for any 
microgravity experiments. 

In the study of dendrite morphology, a computation scheme was 
developed which simulates the evolution of a needle-shaped 
crystal in an undercooled melt in the presence of convective heat 
transfer. The algorithm provides a model of crystal growth which 
includes the physical processes currently thought to govern 
growth. The only approximations are those necessary to the 
numerical solution of the equations, i.e., representation of the 
solution at a finite number of "points" using a boundary integral 
method for tracking the interface. We did not resort to 
linearization nor did we assume quasi-static behavior of the 
temperature field. Thus, the evolution of the interface can be 
tracked in space and time to ascertain its form and stability. 

The algorithm was used to study the influence of convection and 
interfacial energy (the Gi bbs-Thompson effect) on growth 
processes. 

A new class of steady state shapes was found for growth in the 
presence of fluid motion. The relation between growth rate, 
undercooling, flow strength, and the other parameters was 
derived. The relation reduces to the well-known Ivantsov form 
when convection is absent. The stability of these shapes, as 
well as those found in the absence of convection, was 
investigated by following the non-linear evolution of the 
interface after it was perturbed. These interfaces were always 
unstable and the well-known tip-splitting instability appeared. 
Adding the effect of interface curvature on temperature (the 
Gi bbs-Thompson effect) produced new interface configurations, 
which were almost paraboloidal. These shapes were always stable. 
These results are contrary to those found in connection with 
either the theory of marginal stability or microscopic 
solvability. The reasons for this are unknown. The algorithm is a 
solution of the full non-linear problem so some di screpanci es are 
to be expected. However, the qualitative difference in behavior 
deserves further study and the numerical algorithm should be 
carefully checked. If the numerical algorithm is indeed error 
free, then currently accepted theories will need revision. One 
paper derived from this part of the study has been published in 
Physical Review A, others are in preparation. 



The work on protein crystal growth was not carried as -far as 

that on dendrite morphology due to the need to develop 
experimental apparatus. In the -first part of this study the 
influence of fluid motion and other transport processes was 
investigated. Theoretical work disclosed that flow processes 
appear to be too weak to slow crystal growth or cause it to 
terminate. Criteria were developed to indicate when diffusion 
rates would begin to influence crystal growth. According to the 
criteria, none of the extant studies on growth kinetics could 
have been limited by transport rates. In each case, observations 
ceased before the crystals had grown large enough for diffusion 
to play a significant role. 

With the results of the theoretical study in hand, we 
proceeded to develop apparatus to measure growth rates of single 
protein crystals. Apparatus was constructed to grow single 
crystals on a sting under carefully controlled conditions and 
record the growth process with a digital imaging system. We 
believe it is imperative to study growth processes quant i tat i vel y 
in systems where crystals are allowed to grow to the size (larger 
than 0. 1mm) where mass transfer effects could be important. These 
s t u d i e s s h o u Id i n c: 1 u d e m e a s u rem e n t s a f s t r u c ture designed t o 
ascertain whether or not there are changes associated with the 
growth rate itself. At present, we are ready to proceed with the 
experimental study and have arranged to collaborate with a 
protein crystal 1 ographer in Princeton ' s Chemistry Department. 


SUMMA RY St CONC LUS ION S 

DENDRITIC CRYSTAL. GROWTH IN THE PRESENCE OF CONVECTION 
Back gr ound 

Extant theories of cr ystal 1 i 2 at i on deal primarily with 
diffusion controlled growth. However, experimental work on the 
crystal 1 i zat i on of a model compound! 13 shows that natural 
convection has a strong influence at low undercoolings. This is 
particularly vexing since low undercoolings are of interest when 
one seeks to establish a correspondence between theory and 
experiment for dendritic growth. Low undercoolings promote slow 
growth and crystals with large tips, which make photographic 
studies easier. Accordingly, it has been suggested that 
experiments in a microgravity environment would enable one to 
test theories under relatively quiescent conditions. Current 
theories omit convective transport , however, and it is clear that 
a general understanding of the role of convection is necessary. 
Indeed, natural convection is always present in the terrestrial 
environment and it might be advantageous to add forced convection 
to alter the growth process. The presence of convective 
transport makes the theory much more complicated, especially when 
the flow is generated by buoyancy driven motions due to the 
inevitable coupling between the equation of motion and the 
temperature (or concentration) field. Nevertheless, a rigorous 
theory which includes convective transport would be extremely 
useful. A theory describing the effects of forced convection is 
a logical step towards developing a comprehensive understanding. 


Summary Of Completed Mork 

We studied situations where forced convection is aligned with 
the crystal axis. The detailed results of that study are 
contained in a PhD thesis by P„ J. Beaghton. A paper on our 
steady state model was published in Physical Review A ( copies of 
the thesis and the paper are appended to this report) and papers 
dealing with the steady state model and our tracking scheme were 
given at meetings of the American Physical Society Division of 
Fluid Dynamics and The American Institute of Chemical Engineers. 
Other papers are in preparation. 

Ivantsov's theory deals with the growth of needle-shaped 
crystals in a pure, subcooled melt but no allowance is made for 
the effects of interfacial energy on the melting point (the 
Gi bbs-Thompson effect) or convective heat transfer. 

Nevertheless, this theory depicts many features of the growth 
process correctly, so it is the logical vantage point from which 
to consider other factors. Recently, the influence of 
interfacial energy ("surface tension") and anisotropy have been 
studied and a theory known as "microscopic solvability" developed 
12 - 63 . 



To investigate the influence of convection, we first looked 
into the steady state growth of an ax i symmetr i c crystal* Using 
the integral equation describing growth of an isothermal crystal 
in the presence of forced convection aligned with the crystal 
axis, we uncovered a new set of steady, sel f -si mi 1 ar solutions 
analogous to those? of Ivantsov, viz., paraboloids of revolution. 
Here the flow is represented by an exact solution to the 
Navi er— Stokes equations in the Oseen approximation; the 
convective terms in the energy equation are taken into account 
rigorously* A relation between the growth rate, undercooling, 
and strength of the imposed flow was derived* It reduces to the 
Ivantsov result in the? absence of convection* We calculated the 
effects of the flow strength on the Peclet number-undercooling 
relation and, as expected, there is a strong effect* To 
establish the relation between growth rate, tip radius, 
undercooling, and flow strength, however, another expression is 
needed, just as was the case with the Ivantsov theory. 

To develop the second relation needed to set the tip speed we 
could have reworked the microscopic solvability theory (assuming 
it is the correct approach) in a form which includes convection. 
We did not do this. First, it is not yet proven that microscopic 
solvability will explain the selection of tip speeds and sizes 
observed experimentally; more experimental results are needed 
with materials having different degrees of anisotropy* In 
addition, it is not obvious that flow alters events on the 
sub mi c: rose: op i c seal e of the cap i 1 1 ary 1 ength i n a si gn i f i cant 
fashion. It may be that microscopic solvability is unaffected. 
This is certainly an important issfie but some understanding of 
the global character of the problem should be obtained first. 
Thus, our efforts are directed towards developing a numerical 
scheme in which convective transport and non-linear effects are 
taken into account so as to enable us to track large scale 
interface motion. 

A rath e r g oner a 1 n u m e r i c a 3. a 3. g o r i t h m w a s dev e loped to compute 
the motion of an ax i symmetric crystal . (The algorithm is 
described in Beaghton's PhD thesis.) For example, with the 
algorithm we can start a needle shaped crystal in the steady 
state configuration determined by diffusion and convection, add 
surface tension and track the system as it moves to its new 
steady state shape. This solution gives the shape of the 
interface and the relation between the various parameters in the 
presence of flow and capillary effects. Then the interface can 
be perturbed and the stability of the new state ascertained. 
Perturbations can also be added arbitrarily during the evolution. 

Our numerical algorithm mitigates, to a substantial degree, 
the prohibitively large cost and storage requirements of other 
PDE solvers when applied to this sort of problem. The technique 
is based on the transformation of the transient 

convective-diffusion equation and the equations of motion for the 
fluid and the interface into a boundary integral problem. 
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Accordinqly, th ra major part of the computational effort is the 

evaluation of integrals. Traditional PDE integration schemes 
such as finite differences or finite elements would require 
discretization of an infinitely large domain and the subsequent 
calculation of internal values of the temperature that are 
increasingly less relevant as one moves away from the solid-melt 
boundary. Our method is superior to the af orement i oned 
techniques since the internal points can be distributed in an 
optimal fashion and efficient schemes are available to evaluate 
integrals. We recalculate the perturbed flow and temperature 
fields as the system evolves. Interface shape is calculated 
explicitly at each time step and regridding avoided. 

Thus far we have: 

Established that the computation method is "stable" and 
"converges" as the size of the time step is reduced and the 
number of points on the interface increased. In several tests 
we set the system on steady state solutions derived from either 
the Ivantsov theory or our convective theory. The system was 
then allowed to evolve (in time) without adding any sort of 
perturbati on . The robustness of the algorithm was apparent 
from fact that the system stayed on the original solution. 

Investigated the stability of the. Ivantsov solution (which 
makes no allowance for interfacial energy or convection). 

S i n c e t h i s c a r i f i g u r a t i o n i s k n a w r t t o b e u n s t a b 1 e t o 
infinitesimal perturbat i ons it hardly comes as a suprise to 
find it unstable? to finite amplitude perturbati ons. 

Investigated the stability of interface shapes present when 
there is forced convection but no surface tension. They were 
unstable and small, finite amplitude perturbations grew without 
bound in the situations studied. Figure 1 illustrates the 
results of calculations in the absence of surface tension. 

Note the presence of "tip splitting". 

Tracked the evolution from the Ivantsov solution to a state 
where the Gi bbs-Thompson (interfacial tension or energy) effect 
was present. We also added anisotropy, in part. At present our 
algorithm is restricted to ax i symmetr i c shapes, which precludes 
exact computation of anisotropic phenomena. Therefore we simply 
added a term analogous to that used in 2-dimensional 
calculations, i.e., the surface tension was made to depend on the 
angle between the local normal to the interface and the crystal 
axis. In the calculations done thus far the effect of this sort 
of "anisotropy" was small. 

Investigated the stability of steady state solutions 
representing the effects of convection and surface tension. In 
the cases studied thus far, the system was stable to finite 
amplitude perturbations. 


ORIGINAL PAGE JS 
OF. POOR QUAUTY 



FIGURE 1 - A diagram showing evolution of a tip splitting instability on an 
axisymmetric crystal. Heat transfer is by diffusion and convection, surface 
tension is absent. The diagram was produced by the numerical algorithm that 
solves the time-dependent integral equation for the evolution of the 
crystal interface. 


We have constructed a new steady state theory which accounts 
for the effect of forced convection on dendritic growth. It 
shows that heat transfer by forced convection has a strong 
influence on quantitative aspects of dendritic growth but the 
qualitative features remain unchanged , i.e. the shape remains a 
paraboloid of revolution in the absence of the Gi bbs-Thompson 
effect . 

A new algorithm has been developed to track the non-linear 
evolution of a dendritic interface. The scheme accounts for 
effects of convection and surface tension with ax i symmetr i c 
crystals. At present, we can study effects due to the strength 
of the external flow, the? degree of supercooling, the amount of 
surface tension, and, in a somewhat ad hoc fashion, crystal 
anisotropy. However, only a relatively small number of 
situations have been investigated to date and we have not yet 
looked at regions very close? to the tip of the dendrite in any 
detail. It should be possible to examine questions of 
"solvability" in a dynamic sense by computing the interface shape 
in more detail near the tip. Furthermore, the code can be 
expanded to cover non-ax i symmetric growth by allowing for 
azimuthal variations in the shape of the interface and the 
velocity and temperature fields. This would allow us to account, 
for anisotropy in a rigorous fashion. 

PROTEIN CRYSTALLIZATION 

Background 


The crystallization of proteins is a key step in determining 
their molecular structure from x-ray crystal 1 ography , but the 
precise conditions under which a newly isolated protein will 
crystallize are unknown and must be found by trial and error over 
a wide range of pH, ionic strength, and protein concentration. 
Because proteins often form crystals which are too small or too 
disordered to diffract well, finding crystal 1 i zat i on conditions 
is no guarantee that the resultant crystals will be suitable for 
x-ray analysis. It has been suggested that growth of suitable 
crystals is sometimes the limiting step in obtaining structural 
inf ormationC73 . Thus, if the reasons for such contrary behavior 
were known, it might be possible to optimize growth conditions to 
produce higher quality crystals for structure determi nat i on . 
Convection, in one form or another, has been observed to inhibit 
growth but results are still fragmentary C83. If convection 
turns out to inhibit growth in general or alter crystals in other 
deleterious ways, then experiments in quiescent environments 
could lead to improved crystals. A microgravity environment, 
where convection and sedimentation are reduced compared to the 
terrestrial milieu, might provide a suitable venue. 

An intensive study was begun by scientists at NASA and at 
several universities to see what advantages might accrue. The 


overall effort is coordinated through Dr. Robert Snyder, Chief 
of the Biophysics Branch at MSFC* As the work unfolded,/ it- was 
found that knowledge of the kinetics of protein crystal 1 i 2 at i on 
is meager compared to that for inorganic crystals. For example, 
until recently there was no phase diagram for lysozyme, one of 
the? most widely studied proteins, or any other protein » Hence a 
considerable part of the overall effort is devoted to very basic 
studies. Our study belongs to this class. 

The work at Princeton was done by Mr. Marshal 1 Grant , a PhD 
candidate. Because the field is moving rapidly, we kept in close 
contact with others working on the problem. An extensive 
presentation was made to NASA scientists and their col 1 aborators 
on protein crystallisation in March of 1987. The purpose of the 
meeting was to set out our plans to insure that we had a viable 
approach and would not duplicate work already in progress. The 
meeting was chaired by Dr. Snyder. Others in attendance were Dr. 
R. Naumann , Chief Scientist at the Space Science Laboratory, Dr. 
Charles Bugg , Director of the Center for Macromolecular 
Crystal 1 ogr aphy at the University of Alabama (Birmingham) , Dr. 
Franz Rosenberger , Director of the Materials Research Center at 
the University of A1 abama (Huntsvi 1 1 e) and members of their 
research groups. Our plans were given strong endorsement. More 
recently (August of 1987), Dr. Marc Pusey of MSFC visited our 
laboratory to inspect the experimental set-up 5 Dr. Ray Salemme, a 
protein crystal 1 ogr apher at E. I. DuPont visited us in September 
to present a seminar and discuss our work. 

The study of protein crystal 1 ization is difficult because 
protein molecules are extremely complex* and there are strong 
intramolecular interactions in addition to interactions with 
solvent molecules and other proteins. It is often difficult to 
determine the state of a protein system because the 
physicochemical data (phase diagrams, activity coefficients, 
diffusion coefficients, state of aggregation, etc.) have either 
not been determined or have not been published. There? are 
conflicts between data reported by different workers. For 

example, at pH 4, 20 U C, and 50 mg/ml NaCl , the reported values 
for lysozyme solubility range from 1.7 mg/ml [91 to 4.3 mg/ml 
[10,111. Pusey and Gernert [121 recently found that the 
solubilities of the orthorhombic, and tetragonal forms are quite 
different, although differences between the molecular structures 
in the two forms are minorC131. This lack of data makes it 
almost impassible to draw general conclusions regarding growth 
behavior and its relation to the quality of the resultant 
crystals. In the absence of verifiable relationships between 
system conditions and crystal properties, the protein 
crystal 1 ographer is forced to rely on intuition and repetition to 
obtain suitable crystals. Some proteins, moreover, have not 
yielded satisfactory crystals despite these efforts. 

The small size and inherent disorder of protein crystals are 
two major concerns. A third, related, point is the question of 
conformation changes upon . crystal 1 ization and their effect on 
protein crystal growth. 


In the remaining part of this section we summarize results of 
an extensive survey made to clarify issues connected with the 
role of transport processes- Then we describe the apparatus 
constructed for our work. 


Summary of Completed Work 

Our initial effort focused on understanding why protein 
crystals grow slowly and terminate growth at relatively small 
sizes- Theoretical investigations were made of s (i) effects of 
fluid shear on protein binding; (ii) association of protein 
monomer in the bulk; (iii) salt gradients; (iv) contaminants; and 
(v) interface kinetics and mass transfer- The gist of our 
findings is described next; the full report is reproduced in the 
appendix - 

Based on the experimental observation that crystals grown from 
stirred solutions tend to be smaller than those grown from 
quiescent solutions, crystal lographers have recently sought to 
explain the small size of protein crystals in terms of various 
forms of convection- In particular, buoyancy-driven flows have 
attracted attention as a disruptive mechanism in protein crystal 
growth- We examined several scenarios wherein convective flows 
might interfere with the normal binding of protein molecules to 
the crystal surface and found that the fluid velocities which 
a r i s e from d en sit y d iff er en c es a r e t oo sma 13. to p r oduce the 
proposed effects. Specifically, under "normal conditions" shear 
from natural convection appears insufficient? (i) to denature 
individual protein molecules or strip protein molecules from the 
crystal surface (the bond strength is too large to be affected by 
the relatively weak drag force); or (ii) to impose a preferred 
orientation on protein molecules at: the surface (rotational 
diffusion quickly eliminates any bias due to shear). 

Association of monomer may reduce concentrations significantly 
in t h e h u 1 k and t h i s w o u 1 c! r e d la c e grow t h rates b u t e x t. a n t 
experimental techniques are not able to resolve the issue 
unambiguously- Salt gradients due to salt rejection at the 
crystal -f 1 ui d interface appear too weak to influence diffusion 
rates si gni f i cant 1 y - Contami nants whi ch adsorb protein may 

reduce the protein concentration in the bulk significantly but 
experimental studies, of contaminant effects are lacking. 

As crystals grown in more-or-less quiescent environments get 
larger, the protein flux due to natural convection becomes much 
greater than the diffusive flux- At this point convective mass 
transfer effects may affect crystal growth- According to our 
calculations, previous work on protein crystal growth kinetics 
was done on crystals which were too small for convective effects 
to be significant insofar as mass transfer rates are concerned. 
Simply stated, it was- shown that with the small crystals used in 
the quantitative studies, diffusion was so rapid that the 
interface concentration remained unchanged at the initial bulk 
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value throughout the growth period,. Thus, the extant kinetic 

d a i: a d o n o t i n d i c a t e w h e t h e r ri a t u r a. 1 c a n v e c t i o n p 1 a y s a r o 1 e * 

Nevertheless, experimental work by PuseyCBll shows an 
unambiguous effect of forced convection on the rate of growth of 
lysozyme crystals. We believe it is imperative to study growth 
processes quantitatively in systems where crystals are allowed to 
grow to the size (larger than 0.1mm) where mass transfer effects 
could be manifest. These stud ices should include measurements of 
structure designed to ascertain whether or not there are changes 
associated with the growth rate itself. 

At this point we turned to the design and construction of 
apparatus to study the relation between growth rate, size and 
structure?. The apparatus for this was purchased and assembled 
using funds from a Shell Foundation Institutional Grant to 
Pr i nceton Uni.ver si ty . 

Our intention is to examine? the effect of crystal size? on the 
growth and quality of a single crystal . Published experiments on 
the kinetics of protein crystal growth have? been confined to 
crystals which are too small to exhibit significant size effects. 
The experimental procedure used to date (which involves 
introducing supersaturated protein solution into a sample cell, 
nucleating crystals on the sides of the container, and 
photographing the growing crystals) has severe 1 imitations. 

First, only those crystals which are properly oriented with 
respect to the camera can be measured- Second, the crystals are 
often crowded so that the effects of neighboring crystals are 
significant and only average growth rates can be obtained. 
Finally, the surface of the sample cell probably alters the 
growth rate so that the behavior of an isolated crystal suspended 
:i. n s o 3. u t i o n i s s t i 1 1 u n k n o wn . 

We designed our experimental apparatus (Figures 2 and 3) to 
avoid the limitations listed above. Individual crystals are 
nucleated on a glass fiber “sting", which is then mounted on a 
mi crotr ansi ator to provide "xyz " motion for positioning the 
crystal in the cell; a single rotator orients the crystal about 
the vertical axis. The sample cell is approximately 2 cm on a 
side to reduce wall effects. Once the crystal is in position, 
digitized images of the crystal can be captured at specific 
intervals using a video camera. This provides electronic, time 
sequenced images of the growing crystal. The capture time is 
1/30 s and as many as six separate frames can be stored in the 
computer RAM and a video frame buffer. Writing each frame to 
disk can take about 10 seconds, but this does not present a 
serious limitation because most of our work will be in the size 
range where the time between images will be considerably longer 
than the disk access time. (This results from the combination of 
slow growth rates and relatively low magnification required to 
examine crystals larger than 100 microns in diameter.) 

The Mitutoyo FS50 microscope and objectives are capable of 
magnif ications between 2X and 100X with a minimum working 
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FIGURE 2 - Experimental apparatus for the protein crystal growth experiments. 
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Growth Measurement Assembly 



General Layout 


FIGURE 3 - Photographs of the experimental equipment for protein 
crystal growth. 



distance of 20.5mm. Thus, growth measurements of a crystal 
suspended in the chamber can be made without difficulty over the 
entire size range of interest. The microscope is mounted on a 
rack and pinion track to allow coarse focusing of the crystal 
image, while the microscope's focus adjustment will be used for 
fine focusing of the crystal image. A Hitachi KP-232 MOS camera 
sends the image to the video digitizing board (Matrox PIP-640) 
where it can be analyzed using an IBM PC-AT. An image of a 0.2mm 
lysozyme crystal which was digitized through a 5X objective is 
shown in Figure 4. Mounting the microscope; track on a rotary 
table gives us the option of slowly rotating the microscope about 
the crystal so that measurements of different. faces can be made 
in sequence. This feature is an improvement over the techniques 
of previous workers in that the growth of different faces of the 
same crystal can now be measured. The current design allows 
rotation through a minimum of 200 degrees of arc. 

Temperature control utilizes a Lauda RMS-6 refrigerated 
circulator, which can control fluid temperature to within 0.01 
degree K. The circulator pumps water through the sample cell as 
part of a cooling loop. (See Figure 5 for a schematic of the 
apparatus. ) Reservoirs of protein/buf f er solution and 
precipitant (sal t. ) /buff er solution are immersed in the 
circulator's reservoir to maintain the temperature of the feed 
solutions. A peristaltic cassette pump supplies fresh solution 
to the growth cell in a once— through arrangement. Different 
concentrations of precipitant, can be supplied by mixing salt 
solutions of different concentrations, while the total flow rate 
can be adjusted by changing the number of channels used to pump 
the solutions and by a flow control valve upstream of the growth 
cel 1 . 

The pH and ionic strength of the solution are monitored 
continuously using a computer interfaced to an Orion system which 
performs the actual measurement. At this time, there are no 
plans for active control of pH and ionic strength during the 
experiment. 


The Current State— Of —Affair $ 

A careful study of the field shows that the protein crystals 
studied by previous experimenters were not large enough to 
exhibit mass transfer limitations, should they exist. Apparatus 
has been assembled which will allow us to study growth in a 
controlled environment and record and analyze crystal size using 
a video imaging system. With this apparatus single crystal 
growth rates can be measured over a range of crystal sizes. 
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FIGURE 4 - A digitized image of a lysozyme crystal on the sting mounted in the 
apparatus shown in Figure 3. 





effluent protein solution 
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FIGURE 5 - A schematic of the experimental set-up for protein crystal growth. 



R E C 0 i Y i ii E N D A T I 0 M S 

CRYSTAL GROWTH IN THE PRESENCE OF CONVECTION 

Recommendations -for this portion of the work revolve around 
testing (and improving) the interface tracking algorithm. The 
work would provide a theoretical framework for evaluating 
experimental results on crystals grown in terrestrial or 
mi crogravi ty environments- This would lead to a computational 
algorithm which includes all the physics currently believed to 
influence crystal growth. 

Specific topics for investigation ares 

1 . Exploration of the capabilities of the current algorithm- 

To date we have not covered a wide range of undercoolings, 
convective flow strengths, surface tension parameters, or 
other fluid properties. Thus, a general study of the 
program's capabilities is in order- In the process we would 
investigate the effect of flow on crystal shape and stability 
i n some detai 1 «. 

2*. Expand the program so as to allow for fine scale 
r e s (3 1 u t i o n of t h e r e q i o n n e a r t h e dendrite tip- 

T h i s i s m a i n 1 y a c\ u. e s t i. o n o f i nc r e a sod c o m p u t. e r s t o r a q e a n d 

execution time. First one should establish how much can be 
accomplished on the IBM 3081 mainframe and assess the 
advantages in moving to a supercomputer. The use of other 
systems would also be explored. Once this has been done we 
will be in a position to investigate "microscopic 
solvability" in the presence of convection. 

3« Expand the algorithm to account for non-a;-; i. symmetric 
crystal shapes , 

Here again the limitation appears to be one of storage and 
execution time. We already do azimuthal integrations as part 
of the boundary integral technique so the extension should be 
relatively straight! orward. This will allow us to include 
effects of crystal anisotropy on the surface tension. 


PROTEIN CRYSTAL GROWTH 

Here we present recommendations on how to determine relations 
between the quality of protein crystals, crystal 1 i z at i on 
conditions, and protein properties. The experimental work would 
utilize apparatus described earlier. The program would establish 
certain fundamental aspects of protein crystal growth and forms 
an essential part of a broader effort aimed at assessing the 
advantages and disadvantages of growing protein crystals in a 


- 19 - 


microgravity environment. 

Three related topics should be investigated. These are: 


1 . The Effect of Convective Mass Transfer on Protein Crystal 
® Growth 

Crystal growth rates should be measured as a function of 
crystal size (which is equivalent to time) in order to 
determine if convection has an inhibitory effect on crystal 
growth. In some experiments the growth environment should be 
0 quiescent, in others forced convection will be present. 


2. Studies on Crystal Disorder 

Crystals of different sizes should be examined using x-ray 
analysis to determine if there is a relation between crystal 
disorder and size. In particular, the evolution of crystal 
• disorder should be studied to see if it can be related to 

growth conditions. Theoretical studies of crystal packing 
would indicate the effect of packing defects on solvation 
stabi 1 ization of the protein crystal . 


R o 3. e o f C o n f o r m a t i o n Chan g e o n P r o t e in Crystal G r o w t h 

The crystallization of a polypeptide, which is known to 
undergo large conformation changes upon crystallization, 
should be studied. The degree of similarity between the 
polypeptide's crystal growth behavior and that of a 
globular protein, which is not believed to undergo drastic 
conformation changes, would serve? as an indication of the 
relative importance of molecular properties in determining 
crystal growth mechanisms. 
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ABSTRACT 


The motion of the freezing front between a dendritic crystal and a supercooled 
liquid is studied using an interface evolution equation derived from a boundary 
integral transformation of the transient convective-diffusion equation. 

A new steady-state theory is introduced that incorporates the effects of con- 
vection in dendritic growth. It is shown that in the absence of capillary effects 
the shape of the crystal-melt interface is a paraboloid of revolution, similar to that 
found in situations where diffusion is the sole heat transfer mechanism. A relation 
between the supercooling, the product of the tip velocity and tip radius, and the 
strength of the flow is derived which reduces to the well-known Ivantsov theory in 
the absence of convection. 

A non-linear interface-tracking algorithm is developed and used to study the 
temporal and spatial evolution of the dendritic interface. The important role of 
capillarity and convection on the interface dynamics is established and the response 
of the interface to finite amplitude disturbances is examined for the first time. Tip 
splitting is identified as the dominant destabilization mechanism in the limit of zero 
surface tension. Finite surface tension leads to interface stabilization, irrespective 
of the magnitude and structure of the external perturbations. Finally, convection 
significantly decreases the magnitude of the freezing velocity. 
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CHAPTER 1: INTRODUCTION 


The interrelation between non-equilibrium systems and complex growth forms 
was recognized long ago [1]. Over the past decade there has been renewed interest in 
the study of such systems and, in particular, processes of pattern formation in phys- 
ical, chemical, and biological systems. Unfortunately, understanding phenomena as 
common as solidification or dendritic growth is hampered by the mathematical 
complexity of the problem and the subtle effects of microscopic mechanisms such 
as capillarity and interface attachment kinetics. Dendritic* growth, the formation 
of branched, tree-like structures, frequently appears in systems where an interface 
advances into a riietastable phase such as a supercooled melt or a supersaturated 
solution. Often a needle-shaped tip propagates at a constant speed while sidearms 
appear continuously along the sides. The mathematical problem resembles the clas- 
sical Stefan problem, where the diffusion equation for the temperature (or the solute 
concentration) must be solved with boundary conditions specified on the moving 
interface, the propagation velocity of which is determined by the heat (or solute) 
flux. On the other hand, important differences exist: 

• The complex interface shape leads to mathematical problem that defies ana- 
lytical solution. 

• Natural convection due to thermal (or solutal) gradients destroys the mathe- 
matical simplicity of the diffusion equation necessitating a considerable increase 
in computational effort. 

• from the Greek word 8eu6pou (tree). 


• The curved interface is not isothermal. Local deviations from the melting 
temperature of the crystal arise from the interaction between the interfacial 
energy and the local curvature of the solid-liquid interface (the Gibbs-Thomson 
relation). The capillary term leads to a singular perturbation problem which 
has only recently been identified [2-5]. 

Ivantsov’s theory [ 6 ] is a cornerstone of our understanding of dendritic growth 
but deals only with situations where fluid motion and surface tension are absent. His 
theory describes a family of uniformly propagating “needle-crystals” in the form of 
isothermal paraboloids of revolution (or parabolas in two dimensions), characterized 
by the single relation between the dimensionless supercooling A = ( Tm — T 00 )c p /Z- 
and the Peclet number p = pV/2a. Here the tip radius is p, the tip velocity is 
V and a is the heat diffusivity of the melt. The difference between the melting 
temperature Tm and the bulk temperature Too is scaled with the ratio of the latent 
heat of fusion, L, to the heat capacitj 7 , c p . 

The Ivantsov family of solutions is degenerate: for a given supercooling there 
exist an infinite number of paraboloidal solutions since only the product of the 
tip velocity and tip radius can be determined. In addition, linear stability analysis 
shows that the Ivantsov paraboloids are unstable to infinitesimal perturbations. The 
dominant destabilization mechanism is tip splitting, whereby the initially smooth 
tip splits into an increasing number of unstable fingers. This indeterminacy can 
be removed by introducing the effects of surface tension. According to the Gibbs- 
Thompson thermodynamic equation (derived in Appendix E), the (dimensionless) 
temperature of the solid-melt interface is given by Tr=A — (do/p)fC, where 1C is the 
dimensionless local curvature and do = 7 Tmc p /L 2 is a capillary length proportional 
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to the solid-liquid surface tension (see Appendix E). The first analysis of the effects 
of surface tension used the Ivantsov solution as a basic state and treated the capil- 
lary term in the Gibbs-Thompson boundary condition using a regular perturbation 
around the zero surface tension state [7]. For values of o = do/{pp) greater than 
a critical value, a*, a continuous family of stable, modified Ivantsov dendrites was 
found. 

The marginal stability hypothesis, developed by Langer and Muller-Krumbhaar 
[7], employs the linearized stability result along with speculations about the role 
of the non-linear effects. Every needle-shaped crystal displays sidebranches, i.e. 
dendritic growth, which is thought to be a result of non-linear processes. Tip 
splitting is supposed to be a phenomenon described by the linear theory whereas 
sidebranching is non-linear and therefore outside extant stability theories based on 
small amplitude (i.e. linearized) analyses. Accordingly, if we imagine a paraboloidal 
crystal growing into a supercooled melt, the operating point of the system (the tip 
radius and the tip speed that correspond to a given supercooling) would be set as 
follows. Suppose o > a*, because either the tip is slender or the velocity small. 
Non-linear effects leading to sidebrances broaden the tip and a decreases. If a goes 
below < 7 *, tip splitting occurs and forms more slender tips. This suggests that the 
natural operating point for the dendrite is determined by the marginally stable 
solution, a — a*, which furnishes a second relation between p and V. Langer [8] 
acknowledges that the hypothesis is theoretically incomplete and a fully non-linear 
analysis is needed to test it. 

Recent studies have focused on the subtle effects of surface tension on the math- 
ematical structure of the problem. The curvature operator in the Gibbs-Thompson 


equation contains higher order derivatives of the interface shape, multiplied by the 
very small coefficient ( do/p ). This leads to a singular perturbation problem with 
respect to the dimensionless capillary length (or the surface tension) and the in- 
terface shape cannot be represented as a regular perturbation about the Ivantsov 
steady-state. Since the capillary length d 0 is several orders of magnitude smaller 
than the tip radius p, the steady-state correction to the Ivantsov solution in the 
presence of surface tension is negligibly small. As a result, the marginal stability 
theory assumes that the smoothness of the Ivantsov solution is not affected by the 
inclusion of the capillary terms. However, a different approach is needed because 
the singular nature of the problem dictates the interface dynamics and the singular 
behavior is omitted from marginal stability. The recently developed microscopic 
solvability theory is implemented through numerical solution of the equations with 
the surface tension term taken into account rigorously near the tip.* The solution 
is joined onto the Ivantsov solution far from the tip, i.e., the surface tension correc- 
tion to the shape is assumed to vanish as the crystal thickens (exponentially small 
corrections to the interface shape far from the tip are neglected). The microscopic 
solvability theory proposes the following: 

a) For finite surface tension the steady-state correction to the Ivantsov solution 
is not smooth at the tip [2-5]. The numerical solution of an integral evolution 
equation shows that the resulting interface shapes have a sharp tip. However, 
cusp-like tips are not permitted by the diffusional kinetics that govern dendritic 
growth. In accordance to the microscopic solvability hypothesis, all admissible 
interface shapes should satisfy a mathematical solvability condition that requires a 

* In contrast to the marginal stability analysis, the interface shape is not repre- 
sented as a regular perturbation about the Ivantsov steady-state solution. 


zero slope at t.he tip. The slope at. the tip decays exponentially fast as surface tension 
approaches zero and the zero-slope-condition is only satisfied by the zero-surface- 
tension Ivantsov solution. This explains the failure of the regular perturbation 
analysis; the capillarity-induced correction is exponentially small and thus cannot 
be represented in terms of algebraic powers of the surface tension. 

b) The microscopic solvability analysis also indicates that some finite degree 
of crystalline anisotropy produces dendrites with smooth tips [4,5]. The Gibbs- 
Thompson relation is modified to include the effects of surface tension anisotropy: 

T r = A-(d 0 /p)[l-f($ u e 2 )]lC, 

where 0,- represents the angle between the outward-pointing normal and the crys- 
tallographic axis i. The function f(0), which in general depends on the geometric 
parameters of the crystal-melt interface, represents a measure of the surface tension 
anisotropy (/ = 0 corresponds to isotropic surface tension). For sufficiently large 
amounts of anisotropy, there is always a discrete set of finite-surface-tension solu- 
tions that satisfy the microscopic solvability condition and thus have smooth tips 
[4,5]. 

c) Linear stability analysis shows that only the fastest growing dendrite is stable 
[9]. All the other members of the discrete set of allowed steady-states are unstable. 
As a consequence, the fastest growing mode represents the operating point of the 
dendrite. 

The success of the microscopic solvability theory stems from its ability to deter- 
mine a unique operating point for the dendrite. Nevertheless, a number of important 
fundamental issues are still unresolved: 


• The sidebranch emission seems to be the result of the non-linear evolution of 


finite amplitude noise which cannot be incorporated in a linear stability analy- 
sis. Kessler and Levine argue that the tip instability in isotropic systems must 
be thought of as a non-linear instability [9]. They also state that sidebranching 
must be understood via a non-linear analysis of the amplification of finite noise 
as the disturbance moves away from the tip. 

• In most of the theoretical work to date, the temperature (or solute concen- 
tration) field has been assumed to respond instantaneously to changes in the 
interface shape, which in turn implies that the diffusion equation can be solved 
in the “quasistatic” limit. This cannot be true for large supercoolings A where 
the diffusion length is comparable to the tip radius or the distance between 
sidebranches. 

• There is as yet no evidence that the critical amount of anisotropic surface 
tension required to produce a unique, linearly stable, steady-state can be forced 
to remain constant as the interface evolves in time. 

• Experiments in the low supercooling regime show that the sidebranch growth is 
orientation dependent and measured growth parameters such as the tip radius 
and the freezing velocity deviate from their pure diffusion predictions [10]. Ex- 
isting dendritic growth models clearly cannot be used to study these important 
convective effects. 

One of the reasons for the lack of a transient, convective-diffusion model of den- 
dritic growth, aside from the computational difficulties to be described later, is the 
remarkable success of the quasistatic Ivanstsov theory. This simple theory remains, 


despite its limitations, a valuable starting point since it accurately predicts the re- 
lationship between supercooling and Peclet number for a wide range of moderate to 
large supercoolings. Furthermore, the shape of a real dendritic tip region is unmis- 
takably paraboloidal as in the Ivantsov theory. Sidebranches are encountered only 
after one moves five or ten tip radii away from the nose of the dendrite. 

On the other hand, even under carefully controlled experimental conditions, 
thermal gradients in the melt generate buoyancy-driven flows which disturb the 
diffusion-dominated temperature profiles. Experiments by Glicksman and Huang 
[10] show how convection alters the morphology of the crystal and imposes an 
orientation dependence on the growth parameters. Figures 1 and 2 depict the tip 
velocity and the tip radius they measured with pure succinonitrile dendrites. The 
deviation from the pure diffusion predictions is clear at low supercoolings. Figure 3 
shows the effect of convection on the morphology of the dendrite. The temperature 
field is severely distorted and an orientation-dependent elimination of sidebranches 
occurs. Although the purely diffusive Ivantsov theory can be applied to predict 
the experimentally measured Peclet number for supercoolings higher than 1°C (or, 
in dimensionless form, A > 0.05), there is a strong deviation from the Ivantsov 
solution for smaller undercoolings (Fig. 4). Given these striking experimental 
results, it is worth considering how convection alters the theory. Convection has 
not been treated in theoretical work to date for several reasons*. The basic reason 
is the computational difficulty. The boundary integral formulation of the problem, 

* McFadden and Corriell [11] recently extended the Ivantsov solution to include 
the effect of a flow' field due to a density difference between the two phases. 
Their analysis, however, does not apply to motion driven by other forces and 
will not account for the orientation dependence of the experimentally grown 
crystals or explain the low supercooling deviation from the pure diffusion the- 
ory. 


used almost exclusively in recent years, loses much of its simplicity when fluid flow 
is added since domain integrals are required to describe the effects of convection 
in the melt. In addition, the shape of the tip region remains paraboloidal, even 
for very small supercoolings, and this may have prompted investigators to neglect 
any heat transport mechanism that might destroy the paraboloidal symmetry of 
the pure diffusion case. 

The primary objective of this work was to incorporate convection in a tran- 
sient, three-dimensional dendritic growth model based on the convective heat dif- 
fusion equation without any ad hoc approximations. A new steady-state theory 
was developed that shows convection is compatible with the paraboloidal shape of 
the interface. This “convective” steady-state solution was then used as a starting 
point for a non-linear interface tracking algorithm. The numerical scheme solves 
an integral evolution equation for the interface shape and the Oseen hydrodynamic 
equation for the flow velocity simultaneously. Representative results indicate that 
capillarity and convection play important roles in the dynamics of dendritic growth. 

This thesis is structurally divided in two major parts: a steady-state analysis 
of the convective effects is presented in Chapter 2 and a non-linear scheme tracking 
the interface evolution is presented in Chapter 3. The details of the mathematical 
derivations can be found in Appendices A-E. 

In Chapter 2, a complement to Ivantsov’s theory is presented to provide a basis 
for more detailed investigations where convection is present. In natural convection 
the temperature field is coupled to the equation of motion through the buoyancy 
term, which makes the problem all but intractable for most interface shapes unless 
one resorts to numerical methods to solve the partial differential equations. Thus, 


a simpler flow, such as forced convection past a paraboloid of revolution, is used to 
express the salient effects of convection on the solidification front.. 

The velocity field used here is that for flow directed parallel to the axis of a 
paraboloidal crystal. An exact solution to the equations of motion in the Oseen 
approximation is used to represent the flow so that viscous and pressure forces' are 
balanced with a small contribution from inertia. In Sec. 2.1 the integral equation 
that represents the interface shape is analyzed. Capillary effects are ignored, so 
the interface is isothermal. It is shown that uniformly-translating paraboloidal 
solidification fronts are admissible solutions. Then, an expression is given relating 
the supercooling to the Peclet number and the strength of the flow. The application 
of the steady-state analysis to the experimental system used by Glicksman and 
Huang [10] is discussed in Sec. 2.2. 

In Chapter 3, a. non-linear interface tracking scheme based on a boundary 
integral transformation of the transient convective-diffusion equation is introduced. 
With this algorithm all the unresolved issues can be addressed and the validity of 
“steady-state” results, such as the macroscopic solvability hypothesis, can be tested 
with a fully transient, non-linear calculation. 

Boundary integral equations have received increased usage over the last few 
years*. These equations are very useful when evolution of a boundary is to be 
calculated, because, in contrast to treatments using partial differential equations, 
the solid-liquid interface discontinuity is embedded in the singular kernel of the 
integral equation and there is no need to divide the domain in two or more regions. 

* Because of the range of problems, reports in the scientific literature are widely 
scattered. However, a useful treatise on boundary integral methods and their 
applicability to the solution of transport problems has been compiled by Breb- 
bia et al. [12]. 


Recent steady-state models of dendritic growth make use of the integral equivalent 
of the diffusion equation to obtain the crystal-melt interface shape.- The formulation 
presented here, although based on similar principles, actually represents a scheme 
for tracking the interface shape as it evolves in time and space. The presence of the 
convective term, which represents the effects of fluid flow in the melt, increases the 
dimensionality of the integral equation and requires special treatment. The problem 
must be solved numerically, and since the integral kernel is a non-linear function of 
the interface shape, Newton’s iteration method is used to obtain the position of the 
interface at each time step. 

The results presented in this thesis cover the moderate and high supercooling 
regimes and apply only to interface shapes that are single-valued functions of the 
radial coordinate. A summary of conclusions follows: 

(i) The tip is unstable to finite amplitude perturbations if surface tension is ab- 
sent, irrespective of flow strength. Given the results from the linearized theory 
without flow, this is not surprising. But it is important to know that the desta- 
bilization mechanism appears to be tip splitting. Noisy interface perturbations 
seem to invariably focus themselves on the tip region and incipient sidebranches 
are not seen. 

(ii) Surface tension, however small, appears to stabilize the crystal interface for 
any flow strength. A steady-state is always reached, irrespective of the struc- J 
ture or the magnitude of the initial perturbation on the interface shape. The j 
growth velocity of the dendrite decreases with increasing surface tension but 
the difference from the zero surface tension case is small. Convection is found 


to significantly decrease the magnitude of the growth velocity but has no qual- 
itative effect on the stability of the interface. 

(iii) Anisotropy does not appear to have an important effect, qualitatively or quan- 
titatively, on the dynamics of the interface. However, fully-three-dimensional 
(i.e., non-axisymmetric) perturbations will need to be tested before a complete 
understanding of the role of anisotropy can be gained. 

(iv) The operating point of the dendrite is not established from the non-linear 
analysis presented here. However, the results of the analysis call into question 
all those obtained with the linearized theories now extant. Consider, first, the 
marginal stability hypothesis. Given the robust stability of shapes investigated 
here, it appears that the marginal stability hypothesis is wrong. The surface- 
tension adjusted shapes are found to be stable to finite amplitude perturbations. 
Thus, the central feature of marginal stability, a critical value of cr, simply does 
not exist when surface tension is considered in the context of finite amplitude 
effects. Second, microscopic solvability appears incomplete, at best. All the 
shapes investigated in this work are stable when finite amplitude effects were 
considered. Yet one feature of microscopic solvability is that only the fastest 
growing, smooth tip is stable. If this conclusion holds when finer resolution 
is considered at the tip, then the operating point selection mechanism is even 
more subtle than the microscopic solvability would suggest. 

The outline of Chapter 3 is as follows. The formulation of the boundary in- 
tegral equation for the interface shape is presented in Sec. 3.1. a. This evolution 
equation can be used to model growth under the most general conditions and is 
directly extendable to solute-transfer governed dendritic growth. An optimal set 



of numerical discretization points, based on Gaussian quadrature formulae, is in- 
troduced in Sec. 3.1.b. and the boundary integral equation is transformed into a 
set of non-linear algebraic equations. The Oseen equation is also transformed in 
a boundary integral equation which is solved numerically for the velocity around 
the perturbed dendrite (the geometric complexity of the dendritic interface posed 
a challenging computational problem that was successfully handled by the integral 
formulation of the hydrodynamic equation). Sec. 3. 2. a covers the numerical tests 
that sho\v the excellent convergence characteristics of the numerical scheme. Rep- 
resentative results showing the profound effects of capillarity and convection are 
presented in Sec. 3.2.b. The role of anisotropy in the growth dynamics is discussed 
in Sec. 3.2.c. In addition, some ideas for future development are presented, includ- 
ing three-dimensional calculations to test the transient validity of the microscopic 
solvability hypothesis and a finer length scale for the tip region at low supercoolings. 

The core of the mathematical derivations is presented in Appendices A through 
D. The appendices represent an integral part of this thesis and are grouped at the 
end of the manuscript for easy access. Appendix A contains the derivation of the 
integral evolution equation for the interface. The next two appendices highlight the 
difficulties associated with the numerical solution of integral equations. A series 
of variable transformations is employed and the integrals are finally transformed 
into finite sums. The temporal discretization of the boundary integrals is presented 
in Appendix B, while Appendix C describes their spatial discretization. The fluid 
velocity around the perturbed interface is calculated in Appendix D. The hydrody- 
namic equation is transformed to an integral equation which is solved numerically 
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using the boundary element method. The Gibbs-Thompson thermodynamic equa- 
tion is derived in Appendix E. Figures and a listing of the FORTRAN code are 
presented in Appendices F and G. respectively. 
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CHAPTER 2: STEADY-STATE THEORY 


A new steady-state theory is introduced that incorporates the effects 
of convection in dendritic growth. It is shown that in the absence of 
capillary effects the shape of the crystal-melt interface is a paraboloid 
of revolution, similar to that found in situations where diffusion is the 
sole heat transfer mechanism. A relation between the supercooling, the 
product of the tip velocity and tip radius, and the strength of the flow is 
derived which reduces to the well-known Ivantsov theory in the absence of . 
convection. 


2.1. Development of the theory 

Consider the steady uniform propagation of an isothermal solid-liquid interface 
with a constant freezing velocity Vi z , as illustrated in Fig. 5. Densities of solid and 
liquid are assumed to be identical. In a frame of reference traveling with the front 
velocity V, the steady temperature field in the supercooled melt is governed by the 
convective-diffusion equation:* 

a& 2 T + V ttt = v • VT, (2.1) 

Z 

where v is a steady flow field that satisfies the incompressible Navier-Stokes equa- 
tions and the no-slip and mass conservation conditions on the solid-liquid interface. 
Since interface motion does not generate convection unless there is a density dif- 
ference between the solid and liquid phase, v = 0 in the absence of an externally 
imposed flow. 

Capillary effects are neglected and the entire solid is assumed to be at its 
melting temperature Tm , while the bulk liquid phase is supercooled at Too. The 

* The convective term V arises from the motion of the coordinate system. 


heat conservation condition at. the interface z = ((x, y) reduces to 

L 

— Vn- Z = -aVT ■ n, (2.2) 

c p 

where n is the outward unit normal with n z as its 5-component (in the zero-surface- 
tension limit, and in the absence of interior heat sources or sinks, the entire solid is 
isothermal). 

Next, introduce a dimensionless temperature T = (T — T^Cp/L and a di- 
mensionless velocity v = v/f/oo, and scale the lengths with 2a/V. In terms of 
dimensionless variables the equations are: 


V 2 T + 2— = 2A v • VT, 

oz 


VT • n 


*=C(z>y) 


2n 2 , 


T 

*=C(*>y) 


= A 


) 


(2.3) 

(2.4) 

(2.5) 


where A = Uoo/V-, U 0 0 represents the characteristic flow velocity. The interface 
z = ^(x, y) and the solid are now at the temperature A = ( Tm — T 00 )c p /L (the 
dimensionless supercooling), whereas the dimensionless temperature goes to zero as 
2 — ► oo. 

The problem stated above can be cast in terms of integral equations. However, 
since the analysis is length}', it is presented in its entirety in Appendix A. Equation 
(A13) is a reformulation of the problem in terms of the integral evolution equation 
for the interface shape £( x,y ): 

/ oo poo 

dx / dy'2G ss (xr,Xr) 

-oo J — OO 

/ oo poo poo 

dx' dy' dz' 2G aa (x r , x') v(x') ■ VT(x') , (2.6) 

-oo j — oo j C(x' ,y') 


# 
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where x = (a-, y,z) and xr = (x, y, ((x, y)). G aa (x, x') denotes the Green's function 
corresponding to steady heat diffusion due to a point source at x in the reference 
frame moving with velocity V. In the pure diffusion problem A = 0, so the temper- 
ature field is represented by an integral superposition of point heat sources along 
the solid-liquid interface £. The boundary temperature A can then be expressed 
in terms of two-dimensional integrals despite the three-dimensional structure of the 
temperature field. This reduction of the dimensionality of the problem makes the 
integral formulation very convenient when investigating interface motion if convec- 
tion is absent. 

The integral representation of the convective-diffusion equation involves inte- 
gration over the entire fluid domain. (In the limit of pure diffusion, the second term 
on the right-hand-side of Eq. (2.6) vanishes). This increase in the dimensionality of 
the integral equation is due to the lack of appropriate Green’s functions for partial 
differential equations with variable coefficients. Nevertheless, a scheme based on 
the integral formulation has many computational advantages over other methods 
even when there is flow in the melt. 

The integral expression is now used to search for uniformly translating inter- 
face shapes in the presence of fluid flow. Experiments in the “convective” regime 
(low supercoolings) suggest that the tip region remains paraboloidal even when the 
characteristic flow velocity is much larger than the freezing velocity. Thus one needs 
to look for a class of temperature fields (and the corresponding flow fields) that sat- 
isfies the integral equation in cases where the interface is a paraboloid of revolution, 
viz., 

| ( 2 - 7 ) 
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Here the Peclet number, p. can be viewed as the dimensionless radius of curvature 
at the tip of the paraboloid. 

The search for interface shapes is patterned after Pelce and Pomeau [1], who ! ® 

t 

used elementary variable transformations to show that the integral 


r°° 

I dx' 

J — OO 

is independent of x = (x,y,z) if both x and x' represent points on the same 
paraboloidal interface £, i.e. x = (x, y, £(x, y)) and x' = (x', y', £(x', y')). Their 
result was used to demonstrate that the Ivantsov paraboloid is a solution of the 
pure diffusion equation and it can also be used to show that the integral 


L 


dy' 2G 3S (x,x') 


/ oo pa 

dx' / 
-oo J —c 


dy' 2 G a8 {x,x') 


i 


is independent of x when x and x' represent points on two different confocal 
paraboloids, such as those represented by Eq. (2.7) for two different Peclet num- 
bers. Using Pelce and Pomeau’s result, one can produce a class of. temperature 
fields compatible with a paraboloidal interface shape. Since the left hand side of 
Eq. (2.6) is independent of the position vector Xp, one must show that the right- 
hand-side can also be independent of Xr under certain assumptions about the nature 
of the temperature and flow fields. The first integral in the right-hand-side of Eq. 
(2.6) represents the contribution of diffusion and has already been shown by Pelce 
and Pomeau to be independent of the interface position vector xp as long as xp 
represents the same paraboloid as Xp. The second integral in Eq. (2.6), which 
comes from the convective term of the governing equation, must be shown to be 
independent of the position vector xp as xp traverses the paraboloidal interface (. 
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Let 


• 

3 



I 


Ic = /°° dx' f°° dy' f°° dz' 2G ss (x r ,x') v(x') • VT(x') 

J OO J — OO ''C( l/ >y') 


be the convective contribution to the integral expression (2.6). Since paraboloidal 
shapes are under study it is convenient to introduce a paraboloidal coordinate sys- 
tem (diagrammed in Fig. 5): 

z -f \/x 2 + y 2 4- 
w = 

P 

< _ -z + \/x 2 4- y 2 4- z 2 > . 

5 — 

p 

y 

ip = arctan — 

v. x ./ 

In this new coordinate system, the interface defined by Eq. (2.7) is represented by 
the surface w = 1; z — *■ oo is equivalent to w — > oo. Conversely, the surface w = a 
represents a paraboloid confocal to w = 1 with a dimensionless tip radius p' = ap. 
The integral I c becomes 


r2n />oo foo i , i 

I c = dip' ds' l dw' p 3 — - — G ss (xr,x') v(x') • VT(x'). (2.8) 

Jo Jo Jw' = \ * 


If 


P 3 ^4^ V ( x ')VT(x') = A(s»') 


(2.9) 


were independent of s' and <p' , (this will be shown to be true for a certain class of 
flows later) then Eq. (2.8) could be rewritten as 


r°° 2 - r 2n r°° n 2 w' 

h = I dw' — A(w') dip' ds'?—G 33 (x r,x'). (2.10) 

J w f =l Jo Jo * 

Once this is done, the integrals with respect to s' and <p' can be rewritten in terms 
of x' and y 1 using the following coordinate system transformation 


/ dip' f ds'?-£-f(s',ip',w' = w') = ( dx' f dy' g{x',y',z' = ({x',y')) 
Jo Jo ^ J—oo J — oo 


19 


where w' = w' denotes that f(s',ip',w') is calculated along a paraboloidal surface 
of constant w\ and g(s' ,cp' ,w') is the equivalent of in the paraboloidal 

coordinate system. Eq. (2.10) becomes 

roo roo rOO 

I c = / dw' A(w') / dx I dy' G ss (xr, x*) , (2.11) 

JlV 1 = 1 J —oo J — oo 


where 


A ( w> ) = 


= = x',y',z' = 


pw 


1 - 


+ y ' 2 

(pw ') 2 


and 

1 

For every value of w ' , x' and x represent two confocal paraboloids with a ratio of tip 
radii equal to w' . The extended result of Pelce and Pomeau therefore applies directly 
to the integrals over x' and y' in Eq. (2.11), i.e., for a given paraboloidal surface 
w' the integral does not depend on the position vector xr as it moves along that 
surface. A final integration with respect to w' shows that the the convective integral 
in Eq. (2.6) is independent of the interface position vector xr (which represents 
the paraboloidal surface £) as long as the quantity A(s' , y?', w') is a function of w' 
alone. 

To identify a situation where A(s' ,<p' ,w') is independent of s' and , consider 
the uniform propagation of a paraboloidal freezing front in the presence of a flow 
field with a far-field uniform velocity — Uooi z in the direction of the axis of the 


xp = (s. y>, w = 1) = 


P 

x. y. z = - 
,y. 2 





f 






paraboloid. Fig. 5. The fluid velocity v satisfies the no-slip condition on the solid- 
liquid interface (note that in the moving reference frame, v does not include the 
uniform velocity —V\ z which is simply an artifact of the coordinate transformation.) 

Uniform flow past paraboloids of revolution (or parabolas, in two dimensions) 
has been studied by a number of authors and Davis and Werle [2] showed that the 
solution to the Oseen equation 


2A d\ 
Pr dz 


= -VP + V 2 v 


( 2 . 12 ) 


is a uniformly valid approximation* to solutions of the Navier-Stokes equations for 
small Reynolds numbers, Re = pUoo/u=2p\/ Pr, where v is the kinematic viscosity 
and Pr = v/a\s the Prandtl number. Wilkinson [3] derived an analytical expression 
for the Oseen flow velocity v = (u w ,u s ,u v ) past a paraboloid in a uniform stream 
parallel to its axis. In paraboloidal coordinates the velocity is 


v w = 


e -A ... e -Au» 


, / - .E 1 (A U .)-.E 1 (A) 1 

% /ST+7 L AVw£,(A) V £i(A) 


V. = 


> fs E\ (Aw) - Ei (A) 

eT( A) 


y/w + S 


V? = 0, 


(2.13a) 

(2.136) 

(2.13c) 


where E\ is the exponential integral of first order and A = Xp/Pr. On the surface 

w = 1 (z = (), v w = v s = 0, whereas for w — * oo (z — ► oo), v w = —1 and v a = 0. 

The introduction of paraboloidal coordinates yields 

d 2 T d 2 T 
W 8w 2 +S ds 2 

-f (1 + pw - Xpy/w(w + s)v w )t^ + (1 - ps - Apv^(uT+7)u 8 )— = 0 , (2.14) 

* This is true only for three-dimensional flows. In two dimensions the Oseen 
solution for flow past a parabola is not a uniformly valid approximation for 
low Reynolds numbers[12). 
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dr 

dw 


= -P- 


W=1 


lu)=l 


= A, 


( 2 . 15 ) 

( 2 . 16 ) 


for axisymmetric temperature fields. Equation (2.14) can be rewritten as 

d 2 T d 2 T „ s dT „ ,dT p\ A . . 

”5 S? + 5 + (! + + <* - ■ *»> -gj = 


(2.17) 


where 


t 


A- 2^7T7] + V^*vf 


^7 • (218) 


For a paraboloidal solidification front A(u;,s) and T(w,s) must be independent of 
s. It is readily shown that T = T(w ) will satisfy Eqs. (2.15-2.17) with the velocity 
field represented by Eqs. (2.13), so here A = A(w). 

Now Eq. (2.17) can be integrated analytically to give the derivative of the 
temperature with respect to the normal coordinate w , i.e., 


dT_ 

dw 


( 

Pre~ A 1 

<p(l -w)(l + A) + 

1 

1 — i 

4- 

t*J ' 
1 

S 
1 


In w 


Pr 


+ [— 2?i(A) + E,(Aw) + E 2 ( A) - Jfc(A«0] j, (2.19) 

where E\ and E 2 are the exponential integrals of first and second order, respectively. 
One more integration gives the surface temperature of the crystal, viz. A and the 
Peclet number p : 


f°° dT f°° [ 

= — I dw — — = / dwpexp{ p(l — tu)(l + A) + 

Jw= 1 & w J U >=1 ( 


-1 + 


Pr 


,-Al 


EM) 


lntu 


Pr 


+ -F77T [~ E M) + EM W ) + EM) - E 2 (Aw)} (2.20) 

-Cuv/iJ 


Recall that A = A p/Pr. 
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At this point a relation has been derived between the supercooling and Peclet 
number for a paraboloidal crystal formed by freezing of a supercooled melt in the 
presence of convection, viz. 

A = A(p,A,Pr). (2.21) 

Two new parameters are involved: a Prandtl number, v/a , and the ratio of the 
velocity of the flow to the freezing velocity, A. The relation is more complicated 
than that derived by Ivantsov for pure diffusion in that the convective velocity and 
the viscosity of the melt are involved. Furthermore, the parameter A depends on 
the supercooling through V. 



2.2. Discussion of results 


The introduction of convective heat transfer into the equations governing den- 
dritic growth has been shown to leave the paraboloidal shape intact when the flow 
structure has a certain form. This form derives from a solution to the viscous 
flow equations describing the interaction of a single axisymmetric dendrite with a 
uniform flow when the flow is slow i.e., in the Oseen approximation. This sort of 
convection leads to a new family of needle-crystals whose growth velocity is appro- 
priately modified. 

To illustrate the degree to which convection alters the relation between Peclet 
number and supercooling, some representative calculations using the properties of 
succinonitrile are presented in Fig. 6. As the figure indicates, forced convection 

increases the solidification rate substantially when the characteristic flow velocity ! 

I 

is large compared to the solidification velocity. For example, at a dimensionless j 

I 

l 

supercooling of 0.002, the Peclet number with a velocity ratio of 50 is almost twice . 
the Ivantsov value. To emphasize that this is a very weak flow by ordinary standards 
we cite some results from experiments at low supercoolings. 

At a dimensionless supercooling of 0.002, Glicksman and Huang [4] found the 
growth velocity to be roughly 0.8 microns per second and a velocity fifty times this is 
only 40 microns per second which, as the following scale analysis shows, can easily 
arise from buoyancy. The actual structure of the dendritic mass that generated 
flow in the experiments is not known but, as noted by Glicksman and Huang [4], 
it is larger than the radius of an individual tip. Accordingly, each dendritic arm 
is immersed in a flow field configured by the entire dendritic mass. If we assume 


that the characteristic length scale for the crystal mass is /, then a representative 
velocity in a weak flow generated by natural convection is 


/ 




Uoo = C 1 Gr y. (2.22) 


Here C\ is a constant and v is the kinematic viscosity. The Grashof number is 


g£ALl^ 


(2.23) 


where g is the gravitational acceleration, /? is the thermal expansivity of the melt, 
and AL/cp is the supercooling. The velocity U 0 © is then proportional to the char- 
acteristic temperature difference and the square of the characteristic length: 

tfoo«— A' 2 , (2.24) 

i/ Cp 


where the uncertainty of the constant C\ has been absorbed by the “effective” 

convection length l. As a result, the absolute strength of the flow field decreases 

j 

with decreasing supercooling. Using the properties of succinonitrile, we find 


U 0 0 = 740 A/ 2 cm/s 

with l measured in centimeters. Accordingly a dendritic mass with a characteristic 
length of a little over half a millimeter would generate a 40 micron per second flow 
at an undercooling of 0.002. Nevertheless, the obvious differences between free and 
forced convection are enough to deter us from delving further into the experimental 
results until the detailed structure of free convection for this situation has been 
worked out. 
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r 


It should be noted that our methodology differs from that of Glicksman and 
Huang [4], who used an alternative expression for the velocity Uoo , derived from a 
mass-transfer boundary-layer analysis: 


U o 


VG~Vj. 


Here the velocity is proportional to the square root of the Grashof number and can 
be written in terms of the length scale /: 


U a 


gpLi 


A*. 


This expression predicts unreasonably large velocities in the case of the the succi- 
nonitrile experiments (of the order of cm/s) and is clearly inappropriate for very 
low Grashof number flows. 

The theory advanced here makes no allowance for the effects of surface tension 
which is known to have a profound effect on the details of dendritic growth. How- 
ever, this theory provides the necessary starting point for the adaptation of more 
detailed theories to include convective effects which, according to the experimental 
results on succinonitrile, are quite important. The new “convective” steady-state 
solution represents the asymptotic form of the corresponding time-dependent so- 
lution at large distances from the tip. Furthermore, capillary effects are localized 
around the tip and thus the asymptotic solution is also valid for the general case of 
non-zero surface tension. As will be described in Chapter 3, the asymptotic solution 
facilitates the computation of the boundary integrals since the integrands reduce to 
their known asymptotic form far from the tip. In the next chapter, we present a 
non-linear scheme which tracks the evolution of the dendritic interface by solving 
the transient counterpart of Eq. (2.6) numerically. 
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CHAPTER 3: INTERFACE TRACKING 


A non-linear interface tracking algorithm, based on the boundary in- 
tegral equation derived in Appendix A, is used to study the temporal and 
spatial evolution of the dendritic interface. The important role of cap- 
illarity and convection on the interface dynamics is established and the 
response of the interface to finite amplitude disturbances is examined for 
the first time. 

3.1. FORMULATION AND NUMERICAL IMPLEMENTATION 

3.1. a. Theoretical formulation 

In this chapter we study the temporal and spatial departure of the freezing 
interface from the steady-state solutions described in the previous chapter. It is thus 
advantageous to solve the convective-diffusion equation in a reference frame moving 
with the velocity Vi z that corresponds to the underlying steady-state. With this 
coordinate transformation, steadily propagating fronts are replaced by stationary 
interfaces and the transient interface motion is studied separately. 

We consider the evolution of a solid-liquid interface T solely controlled by the 
heat (solute) transport in the solid and liquid phases. The temperature (solute 
concentration) field T obeys the transient convective-diffusion equation, written in 
a coordinate system that moves with a constant velocity V\ z parallel to the z— axis: 

|T + v.V x j-aVjt— v|L = 0, (3.1) 

/ 

where the x vector is represented by the triad (x, y, z)=(r, z ) in cartesian coordinates 
and (r, </?, z) in polar cylindrical coordinates, a(x, t) is the local thermal diffusivity 
of the medium, and v(x, t) is the local velocity in the liquid phase. 
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According to the two-sided “symmetric” model, proposed by Langer [1], ther- 
mophysical properties of the two phases, such as the thermal diffusivity and the 
density, are equal (i.e., ai = as and pi — ps) and independent of temperature. 
Although this is not necessarily true for most materials, succinonitrile, a “plastic” 
organic crystal used extensively for dendritic growth experiments [2], possesses al- 
most identical thermophysical constants in both liquid and solid phases. We adopt 
the two-sided model assumption throughout this chapter but have also derived a 
similar formulation for the one-sided model*. The motion of the interface is related 
to the rate at which the latent heat of solidification is removed away from the inter- 
face through the heat conservation condition at the interface. The latter relates the 
“jump” of the local temperature gradient to the normal component of the interface 
freezing velocity: 

A ■ ^xT(]iq U id) — n • V x T (solid) = ~ c a V ~f~ ~q^ (3*2) 

where the. first term in the right-hand-side arises from the motion of the coordinate 
system and the second, time-varying, term represents the freezing velocity of the 
interface z = ^(r, t) with respect to the moving frame of reference. Here L is the 
latent heat of fusion per unit volume and c p is the specific heat (the symmetric 
model hypothesis implies that the specific heats of the two phases are identical). 

The interface temperature is determined thermodynamically by the Gibbs- 
Thomson relation (see Appendix E) 

T t = T m l_A*[l- e9 (»i, «,)]£{((?,<)}], (3.3) 

* The one-sided model, briefly discussed in Appendix A, assumes no heat flow 
in the solid (the “convective” steady-state solution, derived in Chapter 2, is 
compatible with both models since the solid is isothermal in the zero-surface- 
tension limit). 
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where Tm is the melting point, d 0 =~fTi\ic p / L 2 is the capillary length (proportional 
to the solid-liquid surface tension 7 ), e represents a measure of the anisotropy of 
the surface tension [3,4], g(S\, 62) is an order unity function of the crystal geometry, 
0 , is the angle between the interface normal n and the crystal axis *, and K. is the 
local curvature of the interface, considered positive when the solid bulges into the 
liquid. The Gibbs-Thomson relation is derived under the assumption of local equi- 
librium (see Appendix E), but it serves as a good approximation under conditions 
encountered in dendritic growth. 

Finally, the temperature field asymptotes to constant values as z — ► ± 00 : 

lim T(x, t) = T-too. (3.4) 


We now scale the lengths with 2a/V and times with 4a/V r2 . Eqs. (3. 1-3.4) then 
become 

^- + 2AvV x T-V 2 T-2^ = 0, (3.5) 

n • V x T (liquid) - n • V x T (solid) = -[2 + C(r,i)]n • i 2 , (3.6) 

T r = A-uK{«r,i)}, (3.7) 

lim T(x, t) = 0 and lim T(x, t) = constant, (3-8) 

Z — ►OO Z — > — oo 

where the temperature is now T = (T —Too)c p /L and the velocity is v = v/^ 00 , ^ = 
Uoo/V is the ratio of the characteristic flow velocity scale Uoo to the velocity V of 
the coordinate system, u = do[\ — ^(^ 1 ,^ 2 )] V/ 2a is the (anisotropic) dimensionless 
capillary length function, A = (Tm — Too)c p /L is the dimensionless supercooling, 
z = C(r,t) represents the interface T, and £(r, t)n • i 2 = (d£/dt)h • i 2 is the normal 
component of the interface velocity relative to the moving reference frame. 



In the melt, the velocity v obeys the Navier-Stokes equations for viscous, in- 
compressible flow and satisfies the no-slip condition on the solid-liquid interface*. 
As mentioned in Chapter 2, the solidification of the interface alone does not gen- 
erate fluid flow and in the absence of fluid flow the velocity v is identically zero 
everywhere in the domain independently of the coordinate system choice. In di- 
mensionless form we have (the characteristic length, time, and velocity are 2 a/V, 
4a/ V 2 , and U 0 © = AV, respectively) 

1 dv 2A „ _ _ _2 

jv Tt + iv v ' V * v = + V * v ' (3 ' 9) 

where Pr is the Prandtl number and v = 0 on the interface T. For certain materials, 
including succinonitrile, the Prandtl number is sufficiently large and the temporal 
term in Eq. (3.9) can be neglected, implying that the velocity field responds instan- 
taneously to changes in the interface shape. However, this is not true for metals, 
where the Prandtl number is quite small. 

The coefficient of the second term in the left-hand-side of Eq. (3.9), which 
represents a measure of the inertial characteristics of the flow field, is also small, 
since the flow velocities that occur in dendritic solidification are generally small. 
Nevertheless, the inertial term becomes important at large distances away from the 
interface and must be included in order to obtain a uniformly valid representation 
of the flow field. As an alternative to solving the full non-linear Navier-Stokes 
equations, we use the linear approximation first suggested by Oseen [6]. In the 
Oseen approximation, the non-linear term v • V x v is replaced by the linear term 

* The assumption of equal densities implies that there is no flow through the 
interface due to volume change. 
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(— i s • V x v) utilizing the far-field velocity (— i 2 ). This linearized equation is a uni- 
formly valid approximation of the Navier-Stokes equation in which the viscous and 
pressure forces are balanced with a small contribution of inertia. In Appendix D 
we present a boundary integral method for the calculation of the flow velocity v, 
based on the methodology first applied by Youngren and Acrivos [7]. However, the 
complexity of the interface shape makes this problem computationally more difficult 
than the unidirectional streaming flows past smooth objects studied by Youngren 
and Acrivos. 

The interface shape z = £(r, t) appears only in the boundary conditions for 
the partial differential equations (3.5) and (3.9). As Eq. (3.6) suggests, it is not 
necessary to know the temperature at arbitrary positions around the solidifying 
front since only the temperature gradient normal to the interface is required to 
calculate the growth velocity. Consequently, if Eq. (3.5) and the corresponding 
boundary conditions could be transformed to give equations for the temperature 
gradient at the interface, then the solution of these equations would give all the 
information necessary from the point of view of solidification. Now it turns out that 
such equations can indeed be derived using Green’s function techniques. In general, 
a non-local* integral equation is obtained which relates the normal derivative of 
the temperature field at the interface to the temperature field in the interior of the 
domain. 

Before transforming the transient convective-diffusion Eq. (3.5) into an integral 
equation, it may be useful to explain the physical background of such equations. The 

* A “non-local” equation for a particular point in space contains not only the 
local value of the field variable (or its derivatives) but also the value of the 
field variable (or its derivatives) in other points in the domain. 


surface of discontinuity (in our case, the crystal-melt interface) acts as a distribution 
of “point sources” and “point dipoles.” In the quasistatic case, and in the absence 
of sources or sinks in the domain, it can be shown that the surface source density 
equals the “jump” between the inner and outer normal gradient of the field (i.e., 
the temperature field) and the surface dipole density equals the “jump” between 
the outer and inner field value. In solidification, the temperature is continuous 
upon crossing the crystal-melt interface implying that no dipoles are needed and 
the interface can be represented solely in terms of point heat sources with a surface 
density determined by Eq. (3.6). Finally, it is possible to incorporate the time 
dependence of the heat diffusion equation into the integral equation. In this case, 
temporal integration introduces the effects of the interface history in addition to 
the non-local effects due to the spatial integration. 

The transformation of a partial differential equation into an integral equation 
requires the use of the corresponding Green’s function. This function represents 
the response of the field variable (i.e., the temperature) to an instantaneous point 
source*. The usefulness of the Green’s function lies in the fact that it is a particular 
solution of the adjoint of the differential equation, and thus it already contains some 
of the characteristics of the desired solution.** A shortcoming of this approach 
is that there are no Green’s functions, at least in a useful form, for differential 
equations with variable coefficients. A prime example is the convective-diffusion 
equation (3.5), where, due to the presence of the convective term v • V X T, no 
Green’s function is known for cases where the velocity varies with position. 

* An instantaneous point source is represented by S(x — x')6(t — t'). 

** Note that Eq. (A2), which is the adjoint of the heat diffusion equation, is 
defined in the inverse time space, represented by the time variable t' . 


We now introduce the Green’s function G(x,x'-,t — t') that represents the so- 
lution to the adjoint of the heat diffusion equation in a moving coordinate system 
and in an unbounded domain (see Appendix A), viz., 


G(x, x'\t — t') = ^- 3 - exp 

[47 r(t — t ')] 2 


r — rf-|-[z — + 
4(< - f) 


, (3.10) 


where H(t — t') is the Heaviside function and G(x, x';< — t') = 0 for t < t'. The 
fundamental solution G(x, x'; t — t') represents the transient spatial response of the 
temperature field to a point source at t = t' and x = x' in the moving coordinate 
system. It is singular at (x = x', t = £'), but continuous elsewhere. 

The integral equation corresponding to Eq. (3.5), the transient convective- 
diffusion equation, is 

T(x,f) = J dt' J dr' |2 + C(r',f') G(x,xr{t')]t — t') 

— 2\f dt' f dr' f dz'G(x, x'; t — t') v(x', t') ■ V x <T(x', t'), (3.11) 

J—oo J 


according to the derivation given in Appendix A. The interface xp is at z = C( r )f) 
and Jd r = / 0 °° r dr f* n dtp. It is worth noting that the first term in the right-hand- 
side of Eq. (3.11) represents a point source distribution along the interface with 
density equal to the discontinuous “jump” of the temperature gradient normal to 
the interface, which in turn is directly related to the interface velocity through the 
heat conservation condition Eq. (3.6). The second integral term, which acts as an 
effective volume distribution of sources, represents the effects of convection. Finally, 
the time integration stems from the transient nature of Eq. (3.5) and represents 
the history of the temperature field. 


We now apply Eq. (3.11) to points on the interface x = Xr = (r.z = £(r, t)) to 


yield 


T(x r ,i) = A - ufC{((r,t)} = J dt' J dr' 2 + <(r\ *')] G(x r (t), x r (t'); t - t') 

— 2A / dt' f dr' / d.z'G(xr(t), x'; t — t') v(x',t') • V x >T(x', t'), (3.12) 

J - oo V 

where Eq. (3.7) has been substituted for the interface temperature T(xp,t). 

Equation (3.12) represents an integral evolution equation for the interface shape 
C(r, <), given the temperature field. Without flow (A = 0), the “convective” integral 
in the right-hand-side vanishes and Eq. (3.12) can be solved to determine the only 
unknown, In the presence of flow, however, the convective integral represents the 
effect of a distribution of point sources in the fluid domain with (unknown) volume 
density v(x, t) -V x T(x,t). Thus we need to evaluate the temperature gradient and 
the flow velocity in the interior of the fluid domain (since v = 0 in the solid.) Taking 
the gradient of Eq. (3.11) for x ^ Xr gives 

V x T(x, t) = f~ dt' f dr' [2 + C(r\ t')] V x G(x, x r (0; * - 0 

J — OO J L 

-2A /*" dt' f dr' r dz'V x G(x, x'; t - t') v(x', t') • V x -T(x', t'). (3.13) 

Eq. (3.13) can be directly used to evaluate interior values of the temperature 
gradient for a given interface shape ((r,t). The evaluation of the velocity field is 
discussed in the next section. 


3.1.b. Numerical implementation 

Now we apply the integral equation (3.12) and the auxilliary field equations 

(3.13) and (3.9) to study the temporal and spatial evolution of a crystal interface*. 

* Our methodology actually applies to any advancing front that obeys similar 
equations. 


As mentioned in the introduction, most investigations of dendritic growth involve 
numerical computations of the shapes of steadily advancing fronts and linear stabil- 
ity. Important “steady-state” theories, such as the microscopic solvability hypothe- 
sis, have not as yet been tested by a transient, non-linear, tracking scheme in which 
the position of the interface is updated continuously. In addition, the temporal 
derivative in Eq. (3.5) has been omitted from all previous work on the grounds 
that its effect is small at small supercoolings A. While this is true, the recent 
work on the microscopic solvability hypothesis has been carried out at moderate to 
high supercoolings, where the temperature field does not respond instantaneously 
to changes in the interface shape. The elimination of this “quasistatic” approxima- 
tion changes the solution method significantly, since now 7 the fully-transient integral 
equation is of the Volterra type in time. While the Fredholm type integral equation* 
that corresponds to the quasistatic case has been studied extensively and is suit- 
able for an eigenmode stability analysis, the Volterra type equation is not amenable 
to such a stability analysis. Furthermore, the non-linearity of the integral kernel 
precludes the use of Laplace transforms that usually offer a useful but cumbersome 
alternative to the normal mode analysis. Finally, the incorporation of convection 
in the integral evolution equation increases its dimensionality and requires an iter- 
ative approach to be explained shortly. Convection has been neglected heretofore, 
although many investigators have noted that hydrodynamic effects are important 
in solidification and experiments have indicated the dramatic effects of convection 
in low supercoolings [2]. 

* An integral equation with fixed integration limits. 


( 


The next step is to develop a numerical procedure, based on Eq. (3.12). that 
corresponds to the solution of the full transient convective-diffusion Eq. (3.5). With 
this, shortcomings of other approaches are avoided and a reliable tool for the study 
of the non-linear evolution of crystal dendrites is set forth. 

From t' = —oo to t' = to , the interface shape, the temperature field, and 
the flow field represent steady-state behavior denoted by the subscript (o). That 
implies that in the laboratory frame of reference the interface is a uniformly moving 
solidification front with a constant freezing velocity equal to that of the moving 
coordinate system and £(r, t) is zero. By setting t = t 0 we can depict the steady- 
state propagation of the crystal-melt interface and examine the effects of convection 
on the existing steady-state theories. This aspect of our numerical work has been 
presented in Chapter 2. On the other hand, a perturbation in the interface velocity 
introduced at t' = to causes deviations from the steady-state; the time integrals in 
Eqs. (3.12-3.13) can then be split in two integral parts. The integrals from t' = — oo 
to t' = to represent the dependence of the current interface shape on the underlying 
steady-state, whereas the integrals from t' = to to the current time represent the 
time history of the perturbed interface. Equations (3.12-3.13) now become 


T(x r ,t) = A - t//C{C(M)} = f dt' f dr'2G(xr(t),xr 0 {t'y,t-t’) 

J — OO J 

+ jf dt'J dr' [2 -I- C(r\ O] G(x r (t), x r (t'); t - t') (3.14) 

-2A f dt' f dr' f dz'G(xr(t),x'\t — t') v 0 (x') -Vx'Tofc') 

J — 00 J */Co(r') 

— 2A f'dt' fdr' f°° dz'G(x r (t),x';t-t') v(x',t') -V x >T(x',t'), 

J to J J C(r',t') 
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and 


1 


V*T(x,i) = J° dt'J dr , 2V x G(x,xr 0 (0;<-0 

+ J " dt'J dr' [2 + C(r # ,0] V x G(x,x r (<');<-0 


(3.15) 


/ to f roo 

dt' / dr' / dz'V x G(x,x';t-t') v 0 (x') -V x -T 0 (x') 

-oo J J (> o(r') 

— 2\f dt' f dr' I dz'V x G{x,x'-,t-t') v(x\*') -V x> T(x',<'), 


where xr 0 represents points on the steady-state interface shape z = Co( r )- Note 
that the interface velocity, £(r, t), relative to the velocity of the advancing steady- 
state front is set to zero in the “steady-state” time interval (— oo,<o)- Also note 
that the Oseen approximation to the Navier-Stokes equation (3.9) (see Appendix 
D) is a “quasistatic” approximation and thus does not involve temporal integra- 
tion. Nevertheless, it is implicitly time-dependent since the boundary shape changes 
continuously with time and, as shown in Appendix D, the flow velocity v(x, t) is 
recalculated at every time step. 

At this point it is instructive to compare Eqs. (3.14-3.15) and their quasistatic, 
purely diffusive counterparts used in previous work [3, 4, 8, 9]. 

• In the quasistatic approximation and in the absence of convection only the 
second integral in the right-hand-side of Eq. (3.14) survives out of the eight 
integrals appearing in Eqs. (3.14-3.15)*. 

• In the steady-state models, the time integration is performed analytically over 
the entire time interval (— oo ,t). In order to see why this is not possible in 
our case, consider the second integral in the right-hand-side of Eq. (3.14). 

• With to replaced by oo. 
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In the steady-state case the interface shape Xr (< ; ) does not change with time 
and is equal to the current interface shape xr(f). This simplifies the analysis 
considerably, since the solution of a Volterra type equation is avoided. In the 
transient case, however, the interface shape is a continuous function of time and 
the Green’s function G becomes a non-linear function of the unknown interface 
shape xr(0- This precludes analytical integration and leads to a non-linear 
Volterra type integral equation and to the associated computational difficulties 
mentioned earlier in this chapter. 

• In the absence of convection, Eq. (3.15) becomes obsolete since no interior 
values of the temperature field are required. This reduces the dimensionality 
of the problem and Eq. (3.14) becomes an boundary integral equation. 

Having enumerated the differences between the steady-state and our “transient” 
problem, it becomes apparent that the latter is clearly harder, primarily because of 
three factors: non-linearity, transient nature, and convection. While the first two 
factors combined preclude the use of traditional stability methods and necessitate 
an interface tracking scheme, convection partially destroys the elegant features of 
the boundary integral approach. The convective term in the diffusion equation (3.5) 
is responsible for the domain integrals in Eqs. (3.14-3.15) that increase by one the 
dimensionality of the integral equations. In the numerical scheme that we present 
in this chapter we have coupled the essential features of the boundary integral 
method with a domain discretization that requires a minimal amount of interior 
points. We were thus able to avoid some of the problems with finite difference 
or finite element methods such as the complicated boundary geometry and the 
long computation times usually required for three-dimensional calculations. The 


“differential” approach has not been successful because of the prohibitively large 
computer requirements for such a transient, three-dimensional convective-diffusion 
problem with a moving interface, since each point in the interior domain carries 
an equal weight. The success of our integral approach is based on the fact that 
the “influence” of the convective temperature field on the interface shape decays 
exponentially fast away from the interface and the corresponding domain integrals 
can be accurately calculated in the vicinity of the interface. Summarizing, our goal 
is to calculate iteratively the interface shape at time t using the evolution equation 
(3.14) that contains the time history and the non-local characteristics of the interface 
and using Eqs. (3.15) and (3.9) to update the temperature and velocity fields. 

The temporal integration in the second and fourth integral terms in Eqs. (3.14- 
3.15) cannot be performed analytically since the interface shape, the temperature 
field, and the flow velocity v are all unknown functions of time. To calculate the 
temporal integrals numerically, the time domain [to,t = t;v) is split in N equal 
intervals. Over each such time interval, [t n _i,t n ), field variables, such as the inter- 
face shape, £, temperature, and velocity, are considered to be constant and equal 
to their respective value at the “nodal” point t E[t n -i,t n )* . The actual location of 
the nodal point within each time interval depends on whether the integral belongs 
to Eq. (3.14) or Eq. (3.15). In the discretization of Eq. (3.14), where the objective 
is to use previously calculated values of all the field variables to calculate the new 
interface shape at time t = t^, we place the nodal points t n at the end of each time 
interval. On the other hand, the calculation of the interior values of the temperature 
gradient and velocity is done iteratively and it requires the simultaneous knowledge 

* The error associated with this approximation can be estimated from the mean 
value theorem. 


of the interface shape. Here nodal points are placed at the beginning of each time 
interval and the requisite information for the calculation of the temperature field 
can be retrieved from prior time steps. 

This temporal discretization of Eqs. (3.14-3.15) now allows us to integrate 
analytically over each time interval. In Appendix B we present the result of the 
time integration of each of the integrals in Eqs. (3.14-3.15), denoted by J* and Hk , 
k = 1,2, 3, 4, respectively. The calculation of each integral requires several steps 
but the final result can be expressed in the following concise manner: 


A — t//C{C(r, f)} = I\ + I2 + h + I4 > 


(3.16) 


V x T(x,<) = Ht+Ht + H. + H,, 


(3.17) 


with the integrals (7i — I 4 ) and (H\ — H4) shown in Eqs. (B5-B12.) 

Each of the expressions Ik (or, if*) involves a two- or three-dimensional spatial 
integral (or, a N term sum thereof,) representing the non-local characteristics of 
the problem at the given time interval. Variables in each integral are denoted by a 
subscript which indicates the time at which the field variables are calculated. This 
is exemplified by the integral 
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that belongs to the N term sum in Eq. (B6). Variables with the subscript ( n ) 
represent quantities previously calculated at time t = t n and are treated as constants 
during the spatial integration. Variables with the subscript (n) correspond to the 
unknown values at the present time, t = tpj, but are otherwise treated in exactly 
the same way. 

Our analysis to this point can be applied to any solidification front that obeys 
the evolution equation (3.14) and the auxiliary equations (3.15) and (3.9). However, - 
in order to derive an optimal spatial discretization, we need to define the basic state 
under investigation. From previous steady-state models, including the convective 
model presented in Chapter 2, it is well-established that there exist solutions to Eq. 
(3.14) that represent uniformly translating needle-crystals with a near-paraboloidal 
shape [3,4,8-10]. Although most of these solutions are unstable [11] or possess cusps 
at the tip of the needle [3,4], they do represent the only known shapes that satisfy 
the steady-state equivalent of Eq. (3.14). Since one of our objectives is to investigate 
the non-linear spatial and temporal evolution of such crystal forms, it follows that 
the asymptotic (spatial) characteristics of this system are determined by the far-field 
behavior of the temperature and flow fields that correspond to these steady-state 
needle crystals. The capillary term in Eq. (3.14) modifies the morphology of the 
interface in the vicinity of the needle tip, but otherwise the long-range behavior of 
the temperature field is determined solely by the heat transport mechanisms, i.e. 
heat diffusion and convection. 

We will use the generalization of the Ivantsov family of steady-state needle 
crystals, derived in Chapter 2, to extract the asymptotic behavior of the integrands 
since, at large distances from the interface, the temperature and flow fields are 


essentiality unaffected by the interface perturbation. This generalized steady-state 
paraboloidal solution, denoted by the subscript (o), will also be used to define a new 
set of spatial variables that map the infinite integration domain onto a unit cube. 

The calculation of the spatial integrals involves integration over the entire r, <p 
plane. In addition, the three dimensional “convective” integrals / 3 , I4, #3, and 
H4 require an extra z integration over the semi-infinite interval [( n (r),oo ). The 
immediate goal is to express these integrals in terms of selected values of the field 
variables in time and space. First, the two-dimensional integrals (ty, I2, H j, and 
H2) are examined by looking again at the representative example borrowed from 
Eq. (B6): 


r*n poo 

I dp' I r' dr' 
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where 


\x N - x n '| = yjr 2 + (r') 2 - 2 rr' cos p' + (Ov(r) - < n (r')] 2 . 


For axisymmetric interface perturbations, the ip' dependence comes solely from 
the cosine term in the expression for |xat — x n '|. It is obvious that the integrand 
becomes singular when r = r' and ip' = 0, since the distance |xw — x n '| between the 
points x/v and x n ' (which appears in the denominator) becomes zero. Nevertheless, 


the integrals exist since this weak singularity (1/ |x;v — x n '|, as |xjv — x„'| — * 0) is 
integrable. To demonstrate this we introduce the new variable 


77 = |xw — x n '| = — ' y[r^YcoTip', 


where 


oc(r,r') = r 2 + (r') 2 + (C;v( r ) - Cn(r'))' 


7 (r, r') = 2rr ; . 


The above integral, which can be written as 
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where h is a non-singular function, becomes 
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The reader should note the similarity between the integral C and inverse elliptic 
functions. In fact, in the special case of a stationary coordinate system, C reduces 
to a non-singular inverse Jacobian elliptic function [ 13 ]. In the general case, the 
integrand in C is bounded by the integrand of the stationary case in the entire inte- 
gration interval and thus its existence is guaranteed. To evaluate the integrals with 
respect to 77, the orthogonal polynomial Gaussian quadrature (based on Chebyshev 
polynomials of the first kind) is applied [ 13 ], e.g., 
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For the t] part of integral C we obtain 


A', 
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where 


y/cii + 7 + x/q — 7 i/g + 7 — -y/a — 7 
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The two-dimensional integrals can now be written as sums of N v one-dimensional 
integrals in r' . 

The interface shapes at time tpj (i.e., Cn) and at time t n (i.e., ( n ), both 
asymptote to the paraboloidal (unperturbed) interface shape £ 0 =(p/2)(l — r 2 /p 2 ) 
as r — ► 00. As a result, |xat — x„'| becomes 


lim |x A t - x n '| = 


fcf 

2 p 


and the asymptotic behavior (with respect to r') of the above integrand is deter- 
mined by the exponentially decaying term 

exp{Cn(r')} = exp 



The new variable 



maps the semi-infinite interval [0, oo) onto the finite interval (0,1], with r — ► oo 
corresponding to q — > 1. 

At this point we are faced with the choice of an optimal discretization for the 
q interval, so that the computation of the integrals is both accurate and flexible. 
The satisfaction of both criteria is not a easy task. Most higher accuracy Gaus- 
sian quadrature formulae require a large number of “weights” and abscissas that 


vary with the number of desired points, whereas equidistant, equal-weight formulae 
such as Simpson’s rule suffer from poor accuracy. A quadrature formula based on 
Chebyshev polynomials of the first kind and of odd order was selected. It not only 
offers a high accuracy, but the weights and abscissas are simple functions of the 
number of points A r r , expressed in analytical form: 

J o f(y) dy = 2TT^ ^2N r + j^~ ( 3 - 19 ) 

where 



The implementation of Eqs. (3.18-3.19) transforms the two-dimensional integral to 
a finite double sum in terms of the values of the integrands at selected points i, k 
in the 5,77 unit square (which is equivalent to the semi-infinite r, ip domain.) 

The three-dimensional “convective” integrals I 3 , I 4 , H$, and H 4 involve an 
additional integration in the z' direction. Again we examine a representative integral 
such as 
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found in Eq. (B8). The discretization of the f dr' integral is identical to that of 
the two-dimensional integrals. For the z' integration the asymptotic behavior of 



the integrand as z' — » oo is used to define a new variable transformation that maps 
the semi-infinite interval [£ n , oo) onto the interval (0, 1]. It is easy to show that the 
integrand, excluding the dot product v n (x') •V x >T ri (x'), decays as 1 / z' for z' — ► oo. 
As it was stated earlier, the behavior of the temperature and the velocity fields 
at large distances from the interface is determined by the corresponding “unper- 
turbed” fields To and Vo. Although the flow velocity reaches a constant value, the 
temperature field decays exponentially as z' — ► oo (see Chapter 2.) The integrand 
thus decays exponentially as exp{— 2(1 + A).?'} for z' — ► oo, where A represents the 
relative strength of the flow field. Define the variable 

£ = exp {-2(1 + A) [z' - Cn(r')]} 

which is equal to unity on the interface z' = Cn( r ') and approaches zero as z' — + oo. 
Next use the quadrature scheme defined by Eq. (3.19) to discretize the transformed 
integral from £ = 0 to £ = 1 at the points £ m , m = 1, 2, . . . , N z , where 



At this point we have a set of variable transformations and quadrature rules that 
transform the spatial integrals into weighted finite sums. Each term in the sums 
corresponds to a point (<7i, »/,-*, £ IT n)*, or 

r, = 

(fik = arccos 
z im ~ C( r «) — 2(1 q. 

* The double indices indicate the dependence of t;,* and £, m on g,-. 
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where a(ri) and /?(rj) were defined in Eq. (3.18). The discretized equivalent of each 
spatial integral is presented in Appendix C. 

The integral equations (3.14) and (3.15) have now been transformed to a set 
of non-linear algebraic equations in terms of the interface shape £ n (rq) and the 
temperature gradient V x T n (r t -, z m ) at selected times t n and locations (r,-, z m ). The 
non-linearity can be traced back to the evaluation of the Green’s function along the 
interface. In contrast to the classical one-dimensional Stefan problem, the (gener- 
ally unknown) curved shape of the interface cannot be represented by a constant 
coordinate surface. This, in turn, results in the explicit appearance of £ (or, xr) as 
the argument of the Green’s function in Eqs. (3.14-3.15). The discretized equations 
(3.14-3.15) can be written in the following concise manner*: 

a - .//qcov.i,,)} = fEE I °w(i( r .-.*«)> 

n=l »=1 k -1 


N Nr 

yy v(r t , Z m , ) • V x T(r,*,Z m ,t n _i)}, (3.20) 
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and the auxiliary equation 
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where the interface is axisymmetric (£ independent of fc), J and W are the non-linear 
algebraic operators in Eqs. (C1-C8), and K, is the curvature operator 

r(" + (C ') 3 + (' 

£{CM)} = - V ^ Lr - (3.22) 

r[l + (Of- 

* The superscripts ( D ) and ( c ) denote contributions from the “diffusive” and 
the “convective” integrals, respectively. 
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The velocity v n _i(r,-, 2 m ) at time t „- 1 is calculated from Eqs. (2.13). (D2). 
and (D12-D13). It represents the contribution of two terms: the “unperturbed’’ 
Oseen velocity, derived in Chapter 2, and a Stokes flow correction, calculated in 
Appendix D and representing the effect of the interface perturbation. Since the 
flow field is a solution of the quasistatic Oseen approximation to the Navier-Stokes 
equations, it satisfies the no-slip condition on the perturbed interface exactly and 
asymptotes to a constant velocity at infinity. The velocity is updated using the 
current interface position and is then used in Eqs. (3.20-3.21) to calculate the next 
interface position. It is worthy of note that in the case of pure diffusion, v = 0. This 
simplifies Eq. (3.20) considerably since the spatial integration is eliminated, while 
Eq. (3.21) becomes obsolete since no interior values of the temperature gradient 
are required. 

Eq. (3.20) represents the interface temperature at r = rj and at time t = tpj. 
By varying the index j from 1 to N r , which corresponds to different radial distances 
from the z- axis, we can obtain N r such equations for a given time t = t jy. If we 
assume that the all the variables at the times t n , n = 1,2 , ,N — 1, are known, 
then the set of N r non-linear algebraic equations 
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can be solved to give the N r unknowns, i.e. the interface shape at the points 
rj, j = 1, 2, . . . , AV (Note that the left-hand-side is expressed solely in terms of 
the unknown interface shape, whereas the right-hand-side contains only previously 
calculated quantities.) 

Finally, by varying the value of N, we obtain a set of equations (3.23) repre- 
senting the history of the interface position up to the corresponding time t n. The 
sequential solution of these equations gives the interface position at t = as a 
function of previously calculated values at t — t„ , n = 1,2, . . . , N — 1 . 

In the next section we will discuss the numerical characteristics of this method- 
ology and present representative results. 



3.2. DISCUSSION OF RESULTS 


The non-linear set of Eqs. (3.23). the auxiliary Eq. (3.21) for the temperature 
gradient and Eqs. (2.13, D2, D12-D13) for the velocity represent the numerical 
approximation of the integral evolution equation for the interface shape. We will 
use these equations to examine the growth dynamics of the solid-liquid interface. 
The importance of quantities like surface tension and convective flow strength has 
already been established in the study of steady-state models. In this section we 
will first show that the interface tracking algorithm based on the above equations is 
numerically stable, convergent, and robust. We will then discuss our results which 
clearly demonstrate the significant role of surface tension and convection in the 
transient evolution of the dendritic interface. 

The primary objective of this work was to develop an exact (i.e. without 
ad hoc approximations) mathematical scheme to represent the combined effects of 
interfacial tension and convective heat transport. This was achieved by transforming 
the governing set of equations into a. tractable numerical scheme. In this section we 
present our current results and discuss several additions and generalizations that will 
extend the utility of our computational scheme. It will be shown that our algorithm 
is an effective tool for the study of spatial and temporal pattern formation. 

The important input parameters are of two kinds. Numerical parameters, such 
as time step size and the number of quadrature points, determine the numerical 
characteristics of the algorithm. Their effect on the solutions is discussed in Sec. 
3. 2. a. The physical parameters (Peclet number, surface tension, flow strength, etc.) 
specify the characteristics of the physical system. The results from our calculations 


at selected points of the parameter space are presented in Sec. 3.2.b. Finally, some 
concluding remarks and suggestions for future work are presented in Sec. 3.2.c. 

3. 2. a. Numerical Characteristics 

This section describes the tests that were carried out to examine the robust- 
ness and accuracy of the numerical scheme. Representative results demonstrate 
the ability of the scheme to track known solutions such as Ivantsov’s solution and 
the new “convective” steady-state solution. The convergence characteristics of the 
algorithm with respect to varying grid sizes and time steps were also studied. Fast 
convergence was achieved with a relatively small number of spatial and temporal 
discretization points. 

The underlying basic state for the interface shape is determined by the Peclet 
number and the strength of the flow field. In the limiting case of no flow (A = 0), 
the basic state is the Ivantsov steady-state paraboloid of revolution. As shown 
in Chapter 2, the paraboloidal shape of the interface is preserved for A ^ 0 if 
the flow field satisfies certain conditions. Thus, the new “convective” steady-state 
solution, derived in Chapter 2, is used as the underlying basic state in the numerical 
calculations (since this solution reduces to Ivantsov’s solution for A = 0, there is no 
need to distinguish between the two). The test cases that were examined cover the 
range of moderate to high Peclet numbers (i.e., p > 0.1). Since the characteristic 
length is 2a /V = p/p , smaller Peclet numbers decrease the resolution of the system 
and require a large number of spatial discretization points. A discussion of future 
improvements in the spatial discretization follows in Sec. 3.2.c. 

The discretization of the integro-differential equations (3.14-15) combines the 
numerical approximation of derivatives and integrals by finite differences and finite 


sums, respective]}'. The time derivative of the interface shape, (, in Eqs. (3.14-15) 
represents the perturbation on the interface velocity and can be approximated by a 
finite difference formula. An equivalent approach is to eliminate the time derivative 
altogether by integrating the corresponding time integrals in Eqs. (3-14-15) by 
parts. The first alternative was chosen since it provides one with the flexibility 
to start with either an initial perturbation on the interface shape or an initial 
perturbation on the interface velocity (the interface velocity is directly proportional 
to the gradient of the temperature field and is thus a good measure of thermal 
fluctuations in the temperature field around the crystal.) 

The spatial integrals in the integral equation (3.14) are discretized using the 
coordinate transformations and the quadrature formulae presented in Section 3.1.b. 
The domain is mapped onto a unit cube and the number of the quadrature points 
in each of the three directions determines the spatial resolution of the system. As 
it was mentioned in the previous chapter, the interface shape, the temperature 
field, and the velocity field that correspond to a perturbed interface asymptote 
to the respective unperturbed quantities at large distances from the tip region. In 
the absence of flow, the Ivantsov paraboloid and the corresponding exponentially 
decaying temperature field represent the interface shape and the temperature field 
at infinite distances from the tip region. When flow is present, the Ivantsov solution 
is replaced by the “convective” steady-state solution presented in Chapter 2. The 
known asymptotic behavior can be incorporated into the integral formulation by 
replacing A in Eq. (3.14) with the steady-state integral expression Eq. (2.6). At 
large distances from the tip the integrands in the right-hand-side of Eq. (3.14) 
reduce to those of Eq. (2.6), since the interface shape approaches its unperturbed 


state. As a result, the two sides of Eq. (3.14) cancel each other far from the 
tip and the integration need only be performed up to the point where there is no 
appreciable difference between the calculated interface shape and the corresponding 
unperturbed state. The minimum number of points required to adequately represent 
the interface shape is identified by examining the point furthest from the tip: if 
the calculated deviation from the unperturbed shape is small (C?(10 -3 )), then for 
r > rjsj r the interface shape can be assumed to be that of the unperturbed steady- 
state paraboloid and no more points are required. 

The first important test of the numerical scheme concerns numerical noise: in 
the absence of an externally imposed perturbation, the numerical solution should 
remain fixed onto the underlying basic state as time evolves. A wide range of flow 
strengths was studied (0 < A < 10) and, in the absence of an externally imposed 
perturbation on the interface shape or velocity, no departure from the basic state 
was observed.* The calculations were carried out to arbitrarily large times, thus 
indicating that the accumulation of round-off error does not generate misleading 
numerical instabilities. This is an important result in itself since the correspond- 
ing basic states are often unstable to infinitesimal external perturbations, including 
random noise. The elimination of numerically driven instabilities enables us to con- 
centrate on destabilizing mechanisms such as thermal fluctuations or interface shape 
perturbations. In addition, the amplification of random noise has been proposed as 
the triggering mechanism for sidebranch growth. Thus, the ability to distinguish 

* Note that A = 1 corresponds to a flow velocity that is equal to the freezing 
velocity of the crystal. For moderate to high supercoolings (or Peclet numbers), 
the freezing velocity of succinonitrile dendrites is of 0(cm/s), while comparable 
flow velocities were not observed [2]. 
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between thermal and numerical noise (the latter being virtually uncontrolled) be- 
comes important in the study of non-linear interaction of noise and sidebranching. 

Several tests were carried out to examine the behavior of the solution as the 
number of points used to represent the interface shape and the field variables were 
changed. When a zero-surface-tension solution was used as an initial condition, 
the computed solution would stay on the initial condition as long as v = 0 and no 
external perturbation was imposed. However, once v is changed, the system evolved 
to a new state to account for the effects of capillarity. The following procedure was 
used to test the convergence of the algorithm: The system was started with the 
zero-surface-tension solution in place, i.e. v = 0. Then v was changed to a non-zero 
value and the system was allowed to evolve to a new state. Figures 7 and 8 show 
the behavior of the tip velocity for different numbers of discretization points in the 
r— and ^-directions and for p = 1. The perturbation velocity at other points on 
the interface has a similar behavior. The results shown correspond to A = 0 and 
v = 0.001. Similar results were obtained for A = 0.1 and A = 1.0, indicating that 
the convergence behavior of the algorithm is essentially independent of the strength 
of the convective flow in this range. 

The point r = 0 is not a nodal point and thus one needs to extrapolate from 
the neighboring nodal points using a local Taylor series expansion for the interface 
shape. The 7 r— point result is fairly accurate and the 12 point result is almost 
indistinguishable from the 30 point case. (The maximum deviation from the 30 
point result is 34% for the 7 point result and only 4% for the 12 point result.) The 
convergence is faster in the (^-direction and here less than 10 quadrature points give 
very accurate results: 23% and 2% maximum deviation from the 12 point result for 



the 4 point and the 7 point case, respectively (see Fig. 8). The above results indicate 
that a small number of quadrature points is sufficient to accurately represent the 
growth velocity both in short and long times and over a variation of several orders 
of magnitude in the velocity. 

While the number of quadrature points (N r ,N v ,N z ) determines the spatial 
resolution of the model, the time step size determines the temporal resolution. 
The use of the time-dependent Green’s function and the resulting time integration 
allows the choice of rather large time steps. As an alternative, one can use a time- 
differencing scheme, whereby the temporal derivative in the transient convective- 
diffusion equation (3.5) is approximated by a finite difference formula. However, 
the partial gain from the elimination of the time integration cannot compensate 
for the substantially smaller time steps required by time-differencing [14]. Figure 
9 shows the tip velocity for different time steps and for A = 0.0. Here again the 
zero surface tension solution was used as the initial condition and the system was 
allowed to evolve to a finite surface tension (v = 0.001) state. The results show 
that the algorithm is rather insensitive to variations of the time step size and, 
thus, large time steps give fairly accurate results. The time step At = 0.1 is quite 
large considering the fact that the system reaches a steady-state at times of order 
unity; however, the maximum deviation of the corresponding tip velocity from the 
At = 0.002 result is only 34%. (Similar results were obtained for A = 0.1 and 
A = 1.0). The larger time steps obviously cannot describe the evolution of the 
interface at very small times. Nevertheless, it is important to note that the final 
steady-state is predicted quite accurately by calculations with time steps varying 
over several orders of magnitude. This allows us to interchange large and small 


time steps within a single run. thus accelerating the approach to the final state or 
increasing the resolution over a particular time interval. 

At this point we reiterate that the global approach to a more-or-less iden- 
tical steady-state, irrespective of time step size or number of quadrature points, 
implies that our solution scheme is both unusually robust and numerically stable. 
Considering the complexity of the mathematical problem and the transient, three- 
dimensional structure of the temperature field, this is a gratifying result. 

3.2.b. Representative Results 

The test cases that will be presented in this section demonstrate the effects of 
convection and capillarity on the evolution of the solid-liquid interface. The gov- 
erning equations constitute an initial/boundary value problem in which the initial 
perturbation on the underlying basic state determines the short-time behavior of 
the system. The numerical approximation does not involve a linearization around 
the basic state and thus there is no limitation on the magnitude of the initial per- 
turbation. At this point the algorithm is restricted to interface shapes z — ( (r,t) 
that are single-valued functions of the radial coordinate r. (In simple terms, we 
are currently “counting” the surface elements that are located at the same radial 
coordinate r only once.) The formulation can be generalized by implementing a 
more general transformation of the surface integral in Eq. (All) to account for 
multiple-valued functions of the radial coordinate. This loss of generality does not 
affect our results as long as the initial perturbations are sufficiently smooth func- 
tions of the radial coordinate. In testing the algorithm and in examining the effects 


of surface tension and convection on the temporal and spatial evolution of the in- 
terface, we have chosen to use the zero-surface tension steady-state solutions as our 
initial conditions. These paraboloidal solutions (Ivantsov or “convective”), derived 
in Chapter 2, represent single-valued, long- wavelength perturbations on the under- 
lying basic states with finite surface tension and can be used effectively in studying 
the role of surface tension. 

The role of surface tension in the dynamics of the interface evolution is very 
important. In the zero surface tension limit and in the absence of flow, any pertur- 
bation, however weak, destabilizes the interface and in very short times tip splitting 
occurs with the tip region splitting into an increasing number of growing fingers. 
Tip splitting is characterized by the superposition of interface perturbations with 
an infinite spectrum of wavelengths. While the short end of the wavelength spec- 
trum can only be represented by multiple-valued functions of the radial position, 
the long end of the spectrum is characterized by smoother, essentially single- valued, 
functions of the radial position. The quantitative validity of the calculations is cur- 
rently restricted to single-valued interface shape functions and thus the growth of 
very short wavelength perturbations must be excluded. However, long-wavelength 
perturbations satisfy the single- valuedness requirement up to sufficiently long times 
and can be used effectively to study the dynamics of the interface. Isolated sur- 
face protrusions, smooth perturbations, and random noise were used to show that 
destabilization via tip splitting occurs invariably in the zero-surface-tension limit. 
Figure 10 demonstrates the effect of an interface perturbation on the stability of 
the interface. The initially smooth perturbation, centered at the tip, destabilized 
the interface into an increasing number of growing fingers. The qualitative features 


of the interface destabilization were found to be independent of the amplitude or 
the structure of the externally imposed perturbation. Tip splitting was identified 
under a wide range of Peclet numbers (from 10 -5 to 10) and flow strengths (from 
A = 0 to A = 1). 

The fast convergence of the algorithm was demonstrated in the previous sec- 
tion. While the interface shape is generally independent of the number of interface 
discretization points for non-zero values of u, it becomes strongly dependent at 
the limit of zero surface tension. This peculiar behavior stems from the special 
characteristics of the tip splitting mechanism and is backed by Langer’s theoreti- 
cal analysis of Ivantsov’s solution [15]. Short wavelength perturbations grow faster 
and eventually dominate the interface shape. As a result, the wavelength of the 
fastest-growing mode is constrained by the spacing between adjacent points on the 
discretized interface. In the theoretical limit of the continuous interface, the fastest 
growing component of the spectrum corresponds to wavelength zero and the result- 
ing interface shape resembles a growing “band”. Although the calculated interface 
shapes are strongly dependent on the spatial discretization, is important to note 
that for a given grid spacing the resulting growth velocities are independent of the 
structure and magnitude of the initial perturbation. 

The most striking result of the computations thus far has been the fact that 
small amounts of isotropic surface tension stabilize the interface. The forced per- 
turbation on the interface shape at t = 0 does not lead to interface destabilization 
as in the zero-surface-tension case. As Figs. 11-12 indicate, finite surface tension 
stabilizes the perturbed interface when no flow is present. The perturbation growth 
velocity ((r, <), plotted vs the radial interface position for different times, finally 


reaches a time-independent, spatially-uniform value, thus indicating that a steady- 
state has been reached. The new steady-state is characterized by a constant interface 
velocity 2 + £ that is smaller than that of the corresponding zero-surface-tension 
state.* 

In the initial stage of the interface evolution (i.e., approximately t < 0.3), the 
perturbation velocity £ increases with time and reaches a maximum value that varies 
along the interface. The stabilizing effect of capillarity is evident in Figs. 11-12. 
For v = 0.020, the perturbation velocity reaches a maximum value of about 20% of 
the corresponding zero-surface-tension value of 2 (see Eqs. (3.2) and (3.6) for the 
scaling of the freezing velocity); The growth of the perturbation is much slower in 
the case of a smaller surface-tension-parameter, i.e. v = 0.001. Figure 12 shows that 
the maximum value of the perturbation velocity £ is about 40 times smaller than 
that of the previous case. The magnitude of the capillary term in the evolution 
equation dictates the degree of interface rearrangement toward the new steady- 
state shape. As surface tension increases, the difference between the initial zero- 
surface-tension and the new finite-surface-tension states also increases and “faster” 
transition is required. This explains the strong dependence of the growth velocity 
on the surface-tension parameter u during the early stage of evolution. 

The zero- surface- tension case is characterized by interface instability via tip 
splitting. Irrespective of structure or magnitude, perturbations seem to focus on 
the tip and unstable fingers replace the paraboloidal tip. In the case of finite 
surface tension, however, tip splitting is absent and the maximum perturbation 

* The underlying basic state for t < 0 is the zero-surface-tension steady-state 
solution, derived in Chapter 2. At t = 0, the surface tension becomes finite; this 
corresponds to imposing an interface perturbation, since the interface shape 
must conform to the new value of the surface-tension parameter. 


velocity occurs at about one tip radius away from the tip (see Fig. 11). Given 
that interface instabilities concentrate on the tip, this result indicates that surface 
tension stabilizes the interface by “dissipating” the disturbance away from the tip. 

As time elapses and the perturbation decays away from the tip, a slower phase 
of evolution sets in. The interface velocity becomes increasingly uniform and the 
velocity maximum moves close to the tip (see Figs. 11-12). At about t = 1 the 
perturbation growth velocity changes sign and finally approaches a constant nega- 
tive value. The growth velocity 2 + £ is thus smaller than that of the corresponding 
zero-surface-tension value of 2. The dependence of the total steady-state velocity 
on surface tension is rather weak (see Figs. 11-12), but the role of capillarity in the 
stabilization of the interface is important. There are two distict effects of surface 
tension on the growth velocity in Figs. 11-12: 

(i) The initial perturbation is a function of surface tension and thus the initial 
stages of evolution must depend on surface tension. This dependence is evident 
for t < 0.3. 

(ii) The stabilizing action of capillarity dictates the transition to the new steady- 
state and determines the value of the correction to the zero-surface-tension 
velocity of 2 as i — ► oo (since y 1 for most materials, the new steady- 
state interface shape and velocity are almost indistinguishable from those of 
the zero-surface- tension case). 

Figure 15 shows that the dimensionless time required for the interface to reach its 
new steady state is about 5, thus indicating that the choice of the characteristic 
time 4a /V 2 in the scaling of the governing equations is appropriate (the final time 
recorded in each of the Figs. 11-14 represents a relative difference from the new 



steady-state growth velocity of less than 10~ 3 ). It is worthy of note that although 
the steady-state velocity differs only slightly from that of the zero-surface-tension 
case, the transition toward this new steady-state is rather complex and cannot be 
captured by a linearized evolution equation. 

Only axisymmetric perturbations have been examined thus far. As a result, 
surface-tension anisotropy depends only on the angle 9 between the interface normal 
and the 2 -axis. We can assume that the anisotropic correction to the surface tension, 
defined in Eq. (3.7), is equal to ei/cos(40), where e represents a measure of the 
anistropy of the material and the factor 4 in the argument of the cosine corresponds 
to the fourfold symmetry of succinonitrile. Figure 15 indicates that in such a case, 
the anisotropic term in the interface evolution equation has a negligible effect on the 
interface dynamics for e as large as unity. However, anisotropy has, by definition, 
a three-dimensional nature and the resulting interface shapes should, in general, 
be non-axisymmetric. A further discussion of anisotropy and three-dimensional 
calculations follows in the next section. 

The combined effects of convection and capillarity on the interface velocity are 
shown in Figs. 13-14. For given supercooling A and flow strength A, the zero- 
surface-tension solution, derived in Chapter 2, defines the underlying basic state for 
t < 0. At t = 0, v changes value and thus a surface-tension dependent perturbation 
is imposed on the interface. The transient response of the interface varies with A, 
but the interface invariably reaches a new steady-state. Figures 13 and 14 show the 
interface velocity as a function of the radial position at different times for A = 0.1 
and A = 1.0, respectively. In comparing Figs. 13-14 with the corresponding Fig. 
12 for A = 0.0 (u = 0.001 and Pe = 1), it is important to note that each figure 


corresponds to a different- supercooling A, given by Eq. (2.20). Since the level of 
supercooling, and not- the Peclet number, is controlled by the experimentalist, Fig. 
6 can be used to obtain the correct A for given Pe and A. However, for A < 1 
the differences are small and a constant Peclet number can be used. Furthermore, 
the flow strength A = Uoo/V depends on the scale factor for the freezing velocity, 
V. Since the actual velocity-tip radius relation is still unknown, V must be kept 
constant if one is to compare interface speeds corresponding to different levels of 
convection. 

As in the A = 0 case, the interface is quickly reshaped to conform to the 
new value of surface tension. The qualitative characteristics of interface evolution 
remain the same, but the coupling of capillarity and convection produces a more 
complex transition toward the new steady-state. Figures 13 (i.e. A = 0.1) and 14 
(i.e. A = 1.0) show a perturbation- velocity maximum of about 1.6% and 13.0% of 
the corresponding zero-surface-tension value of 2, respectively (the A = 0.0 value is 
0.5%). In addition, a comparison of Figs. 12-14 suggests that the “dissipation” of 
the disturbance away from the tip becomes slower as A increases (note the increasing 
number of local velocity maxima in the t = 0.2 curves). This non-linear interaction 
of the flow field and the interface disturbance seems to be linked to the orientational 
dependence of the interface morphology shown in Fig. 3 (the rigorous analysis of this 
phenomenon requires a (^-dependent flow field and is discussed in a later section). 

After the initial reshaping, the interface slowly approaches its new steady-state, 
characterized by a negative correction to the zero-surface-tension interface velocity 
of 2. Figure 15 shows the decay of the perturbation velocity at the tip toward its 
steady-state value. For A = 1.0, ((0,< = 0-3) and ((0,t °°) are equal to 5.5% 


and (7.5%) of the zero-surface-tension steady-state value of 2. For A = 0.1, these 
values are about one order of magnitude smaller. As in the case of A = 0.0, the 
time required for the interface to reach its new steady-state is approximately 5. 

In conclusion, the total steady-state interface velocity decreases monotonically 
with respect to both increasing surface tension and flow strength. However, the cal- 
culations for v = 0 suggest that convection does not contribute to the stabilization 
of the interface. In the following section some conclusions are drawn and several 
ideas for future development are presented. 

3.2.c. Conclusions and future work 
1. Conclusions 

The primary results of this chapter are the following: 

(i) A non-linear tracking algorithm, based on the integral transformation of the 
transient convective-diffusion equation, has been developed for the study of 
complex growth forms. The ability to treat finite amplitude disturbances and 
the applicability to arbitrary geometries make this numerical scheme a valuable 
tool in the understanding of pattern formation. Analogous formulations have 
been used in the past for the study of dendritic growth and viscous fingering 
in Hele-Shaw flow cells [3,4]. The most important differences between this 
analysis and the ones preceeding it are: 

1. The temporal derivative in the convective-diffusion equation has been ne- 
glected in previous quasistatic models. Howevr, scaling arguments show 
that at large supercoolings this term becomes imprortant and cannot be 


ignored, at the expense of time integrals that are quite cumbersome to 
calculate. 

2. The study of the convective effects on the interface evolution requires the 
calculation of spatial integrals over the entire liquid region, absent from 
similar work on diffusion-controlled growth. The algorithm employs a se- 
ries of variable transformations and an optimal set of Gaussian quadrature 
points that evolves in time so that the temperature field is calculated at 
the minimum number of interior points. 

3. The interface shape appears in the non-linear kernel of the integral equa- • 
tion. Previous efforts to linearize the evolution equation and carry out a 
linear stability analysis do not capture the non-linear features of dendritic 
growth. Sidebranching and tip splitting can be thought of as results of the 
non-linear interaction of interface perturbations and capillary forces. The 
algorithm presented here applies the Newton’s iteration method to solve 
the non-linear evolution equation as a function of time and space. It can 
thus be used to address important issues such as the amplification of finite 
amplitude noise and the origination of sidebranches. 

(ii) The numerical characteristics of the algorithm were tested as a function of the 
spatial and the temporal discretization of the interface. Fast convergence was 
achieved with a small number of quadrature points. In addition, the solutions 
are rather insensitive to the time step used for the discretization of the tem- 
poral integrals. Finally, in the absence of an externally imposed perturbation, 
the algorithm tracks the underlying basic state and numerical noise does not 


introduce artificial instabilities. 


(iv) The tip was found unstable to finite amplitude perturbations if the efFects of 
capillarity are ignored. The dominant destabilization mechanism is tip split- 
ting. The tip region degenerates into a increasing number of unstable fingers 
and interface perturbations of arbitrary amplitude seem to invariably focus 
themselves on the tip region. This result, which is independent of the flow 
strength, represents the first evidence that tip splitting is not restricted to 
infinitesimal perturbations. 

(v) Surface tension, however small, appears to stabilize the crystal interface for any 
flow strength by “dissipating” the disturbance away from the tip. A steady- 
state is always reached, irrespective of the structure or the magnitude of the 
initial perturbation on the interface shape. The growth velocity of the den- 
drite decreases with increasing surface tension but the difference from the zero- 
surface-tension case is small. 

(vi) Anisotropy does not appear to have an important effect, qualitatively or quan- 
titatively, on the dynamics of the interface. However, fully-three-dimensional 
(i.e., non-axisymmetric) perturbations will need to be tested before a complete 
understanding of the role of anisotropy can be gained. 

(vii) Convection reduces the total growth velocity of the new steady-state, but does 
not seem to contribute to the stabilization of the interface. In the zero-surface- 
tension case, convection does not affect the qualitative features of tip splitting. 
However, the role of convection in the low supercooling regime needs further 
investigation and is discussed later in this section. 

(viii) The operating point of the dendrite is not established from the non-linear 
analysis presented here. However, the results of the analysis call into question 
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all those obtained with the linearized theories now extant. Consider, first, the 
marginal stability hypothesis. Given the robust stability of shapes investigated 
here, it appears that the marginal stability hypothesis is wrong. The surface- 
tension adjusted shapes are found to be stable to finite amplitude perturbations. 
Thus, the central feature of marginal stability, a critical value of a, simply does 
not exist when surface tension is considered in the context of finite amplitude 
effects. Second, microscopic solvability appears incomplete, at best. All the 
shapes investigated in this work are stable when finite amplitude effects were 
considered. Yet one feature of microscopic solvability is that only the fastest 
growing, smooth tip is stable. If this conclusion holds when finer resolution 
is considered at the tip, then the operating point selection mechanism is even 
more subtle than the microscopic solvability would suggest. 


2. Future Work 

The first extension of this work should be toward the study of three-dimensional 
(i.e. non-axisymmetric) interface perturbations. Since the formulation is applicable 
to fully three-dimensional systems and since the modification to anisotropic surface 
tension is straightforward, this numerical scheme can be used to test the microscopic 
solvability theory under the most general conditions of transient interface growth. 
The validity of the microscopic solvability theory at finite times can also be tested 
in a different way. Instead of prescribing the amount of surface tension and the 
corresponding anisotropy and then examining the slope at the tip, one can do the 
opposite: fix the tip slope to be zero and then calculate the time dependence of 


the capillary parameters. The algorithm requires the solution of a set of non- 
linear equations with an equal number of unknowns. A logical choice for the set of 
unknown variables would be the interface shape function ( calculated at the various 
quadrature points. However, there is no mathematical restriction applying to the 
replacement of one or more of these unknowns with an equivalent number of the 
physical parameters of the system. Therefore, by fixing the slope at the tip, one of 
the unknowns is effectively eliminated and a new degree of freedom is introduced in 
the physical parameter space. By allowing the amount of surface tension anisotropy 
(or, the amount of surface tension itself) to vary, one will be able to monitor the 
microscopic solvability condition through time and thus test its validity in transient 
dendritic growth. 

The accurate calculation of the slope requires a fine discretization of the tip 
area. At the present time a high density of points close to the tip cannot be pre- 
scribed without introducing an unnecessarily large number of discretization points 
in the radial direction. This stems from the nature of the Gaussian quadrature 
formulae and the associated fixed distribution of points within a given interval. Fu- 
ture refinements of the algorithm should include a choice of “specialized” quadrature 
schemes in addition to the all-purpose Gaussian quadrature scheme presented here. 
Such flexibility would also speed up computations and reduce storage requirements 
in the case of an initial perturbation with a narrow range of wavelengths. Calcu- 
lations for long wavelength perturbations require a “larger” domain but make no 
use of the increased resolution of the tip region, with the opposite being true for 
short-wavelength perturbations localized at the tip. 



Experimental evidence indicates that the effects of convection become impor- 
tant at very low supercoolings (A < 0.05) [2]. All the computations, however, were 
done at Peclet numbers greater than 0.1 which belong to the “diffusive” regime 
(see Chapter 2). The study of convection at lower Peclet numbers (or supercooling 
levels) requires a higher spatial resolution around the tip. In Sec. 3.1. a the length 
scale is defined as 2a /V. Using the definition of Peclet number ( p = Vp/ 2a, where 
p is the tip radius of curvature), the above scale becomes p/p. For small Peclet 
numbers, it is obvious that this length scale is much larger than the tip radius 
and the spatial resolution in the tip region is severely decreased. Therefore, at low 
Peclet numbers, where convection is expected to play a stronger role, the tip radius 
p should be used as the more appropriate length scale near the tip. Also, values of 
A greater than unity are possible in the low supercooling regime where the freezing 
velocities are quite small. It thus remains to be seen whether the use of a finer scale 
will show that convection qualitatively affects the stability of the freezing interface 
at low supercoolings. 

Finally, the convective steady-state solution, derived in Chapter 2, applies only 
to axisymmetric flow fields, oriented along the z— axis. The study of the orientation- 
dependent elimination of sidebranches, shown in Fig. 3, requires the use of more 
general (i.e non-axisymmetric) flow fields. The integral formulation is not restricted 
to axisymmetric flow fields and can thus be used to examine growth subject to an 
arbitrary flow field (at the expense of an increase in the dimensionality of the 
equation and in the computation time). 
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APPENDIX A: 


Derivation of the integral equation 


In this appendix we derive the integro-differential equation that describes the 
evolution of the dendritic interface. The transient convective-diffusion equation is 
first written in differential form: 


f + v.V s f-*V!f-V§ = 0, 


(Ala) 


where x = (x,y, z) = (r, z) = (r,tp,z), a is the thermal diffusivity and Vi* is the 
constant velocity of the moving coordinate system. The hydrodynamic velocity field 
v satisfies the no-slip condition on the interface and is zero in the solid. Under the 
assumptions of the two-sided model proposed by Langer [1], the thermal diffusivities 
and the densities of the two phases are considered constant. We scale lengths with 
2a/V and times with 4a/V 2 : 


s) r r F) r r 

Ij- + 2A v . V„T - V’T - 2^. = 0, 


(Alb) 


where A = Uoo/V is the ratio of the characteristic fluid velocity scale to the velocity 
of the coordinate system and T = (T — T^Cp/L. The temperature field vanishes 
as z — ► oo and is also constant as z — * — oo: 


lim T(x, t) — ♦ constant. 

z— *±oo ' ’ 7 


(Ale) 


The temperature also satisfies the thermodynamic boundary condition 


Tr = A-uAC{C(r,t)} 


r 


with A = ( Tm — T^Cp/L and v = doV/2a ; K, is the dimensionless curvature. 
Finally, heat conservation at the interface requires that 


n • V x T (liquid) - n • V x T( soUd) = -[2 + C(r,*)]n- i 2 . 


(Ale) 


We now introduce the fundamental solution G(x, x';t — t') of the transient 
diffusion equation in an unbounded domain; in the moving coordinate system, it is 
defined by 

Bn _ Bn 

(A 2) 


d< 7 -Vl,G- 2 = -S(x - x') S(t - t'). 


dt 


and G(x,x']t — t') = 0 for t < t'. The fundamental solution (or Green’s function ) 


is 


G(x, x ; t — t ) = j- exp 


r-r'| 2 + [z-z’ + 2(t-t')} 


Ml* 


, (A3) 


f 4 7r(t — <')]* * t 4(<-<0 

where H(t — t') is Heaviside’s step function [1]. By taking the Fourier transform of 
Eq. (A2), Langer [1] derived the following integral representation of G(x, x'; t — t'): 


, lx f°° du f dk 
G(x, x ; t — t ) = / — / 

J-oo 2tt J (27 r 


00 du f dk e M<-0+«'Mx-x') 
(27r) 3 io) + — 2 ik z 


(A4) 


Following the notation of Caroli et al. [2], we define: = lim c —o+ ( t — c). Then 


G(x, x'\t — t-) = 6(x — x'). 


(A5) 


Let k = q + fc 2 i 2 . Eq. (A4) then becomes [2] 


G(x,x';i-i') = J°° 

dq e <q < r “ r '> 


(A6) 


I 


(2 tt) 2 2[m(g,w) - 1] 


exp (- {z - z' + \z - z'\{m(q, a/) - 1]}) , 
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where 


m(q, u) = 1 + (1 + iu -f q 2 ) : 


with Re(l + iu> + g 2 ) 2 >0. 

We now multiply Eqs. (Alb) and (A2) with G(x,x';< — t') and T(x',f'), re- 
spectively, add and integrate in time and space: 


Jj* J M £jz‘ {g(x,xm - f)^£l + T(x-,o aG(x ’^; ' - ' ■ > } + 

r~ dt' f dr' r dz' {-G(x,x';« - «')V*,T(x',i') + T(x\t')V*,G(x,x';t- <')} + 

*/ — oo */ */ — oo 


/ d/' / dr' f dz 1 1 -2G(x, x'; t - t') 

J — oo «/ */ — oo l 


0T(x',f') 


f^-2T(x\*')^^ 


02' J 


[*- dt' f dr' r dz'2X G(x, x'; i — <')v(x', f') • V x »T(x', t') =0, (A8) 

J —oo J 

where the lower limit for the 2 ' integration in the last term is £(r',<') since the 
velocity v is zero for 2 ' < C(r 

By grouping similar derivatives together and applying Green’s 2nd theorem to 
the second term in the left-hand-side of Eq. (A8) we get 

J dt' J dr' J dz' ~ [G(x, x'; t — t')T(x', t')] + 

J dt'jdT'ii' • |g(x, xr'{r' ,t')\t — t')V X 'T(x.' ,t') 

- T(x' , t')V X 'G(x, xrV, t'); t - <')] • - 
J dt' J dr' J dz'2-^-j [G(x, x'; t — t')T(x', t')] + 

f ' dt' f dr' f°° d 2 ' 2 AG(x,x';<-t')v(x',<')- V x .T(x',<') =0. (A9) 

J —00 J ^(r', t') 


The second term represents integration over a surface enclosing the interface and 
placed at an infitesimal distance from it. 
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Next the z' integration in the third term of Eq. (A9) is performed so that the 
term becomes 


r dt> i 

1 

f — oo J 

L 


lim G(x;r',zi;< - <')T(r',zi,<')|. 

[ — ► — OO J 


From Eq. (Ale) we know that the temperature as z — ► oo is zero and thus the first 
term in the above expression vanishes. The temperature reaches a constant value 
at z — ► — oo and can be taken out of the integral sign so the remaining integral is 


f dt' [ dr' lim G(x; r', z 1 ; t — t'). 
J- oo J Z ^°° 


Using Eq. (A7) and the fact that lim c _o+ G(x, x'; — c) = 0, the integral can be 
written as 

** / dr' / «•“('-' ) / -ZL -- — — — — ; = o. 

J J- oo 2* J (2ir) 2(m(9,io)-l] 

Thus, the third term in Eq. (A9) does not survive even if the the temperature at 
negative infinity is non-zero. 

We can now perform the time integration in the first term of Eq. (A9): 

f dr' f dz' G(x, x';< — t_)T(x',t_) — lim {G(x,x';t — t')T(x',t')} . 

J J - oo L t'—-oo 

With the help of Eq. (A5) and of the fact that lim^-o-oo G(x, x'; t — t') — 0, this 
integral becomes 

J dx'6(x - x')T(x',<_) = T(x, t), (A10) 

where <_ has been replaced by t since the temperature is a continuous function of 


time. 


Returning to Eq. (A9), we find 


T(x,<) = 

j d< J dr'll ' ■ [G(x,x r '(r',( , );i-(')V lc .r(x',(') 

- T(x', i')V x ,G(x, x r '(r', f')i (-(')] + 

TOO 

dr' dz'2\G(x,x'-,t-t')v(x',t')-V xl T(x',t'). (All) 

A(r'.t') 

Since the domain fi is infinite, we assume that it is bounded by a surface that has 
been removed to infinity. Langer [1] has shown that in the absence of external fluxes, 
the contribution of this outer surface to the second term of Eq. (All) is zero. We 
thus restrict our attention to the crystal-melt interface I\ The surface integral in Eq. 
(All) contains the boundary values of the temperature and its normal derivative 
and the continuity of the temperature implies that the term proportional to the 
gradient of the Green’s function cancels out. We can now evaluate the remaining 
term, proportional to the normal component of the temperature gradient, using Eq. 
(Ale): 



T(x, t) = J'~ dt'J dr' [2 + C(r', *')] G(x, x r '(0; t - t') 

- 2A f dt' f dr' f dz'G(x , x'; t - t') v(x\ t') • V x .T(x', <'), (A12) 
J — 00 J J C( r'.tM 


since for single-valued z = C(i*)0 it is true that 


n • i 2 dT = dr = r dr dip. 

Equation (A 12) is the integral equivalent of the transient convective-diffusion 
equation in the moving coordinate system. In the case of a steadily propagating 


interface, i.e. with a uniform velocity equal to that of the coordinate system, Eq. 
(A12) reduces to the steady-state equation 


T(x) = 2 / dr'C? as (x, xp') — 2A / dr' j dz'G ss (x, x') v(x') • V x <T(x'). (A13) 
J J J C(r') 

Here G 93 (x,x') represents the Green’s function that corresponds to steady diffusion 
in the moving coordinate system and is given by 


G„(x,x') = f dt 1 G(x,x',< - t'). 
J — OC 


(j414) 


The integral equation that corresponds to the one-sided model is identical to 
Eq. (A 12) at the limit of zero-surface- tension since there is no flux through the 
isothermal solid. In the presence of surface tension, the convective effects are rep- 
resented by a term identical to that of Eq. (A 13). On the other hand, the diffusive 
contribution is with the one sided model is different since heat cannot flow' through 
the solid (as = 0). The derivation of the corresponding “diffusion” term has been 
presented by Caroli et al. [2]. 
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APPENDIX B: Temporal discretization 


In this appendix the time discretization of the integral evolution equations 
(3.14-3.15) is described. 

T r = A - i/£{C(r,<)} = / dt' j dr' 2 + C(r',f')j G(x r (t),xr (<');* - t') 

J — OO J 

— 2A f dt' f dr' f dz'G(xr(t),x';t — t') v(x',t') -V' x T(x',t'), ( B1 ) 

J — oo J J C(r',f') 


V*T(x,<) = J dt' J dr' [2 + <(r', *')] V x G(x, x r '(t'); t - t') 

- 2A [*~ dt' f dr' f°° dz'V x G(x , x'; < - <') v(x', t') • V' x T(x\ <')• (52) 

From t = —00 to t = to) the interface shape, the temperature field, and the flow 
field are represented by steady-state solutions denoted by the subscript ( 0 ). At time 
t = to a perturbation on either the interface shape, or the temperature field, or other 
boundary conditions causes deviation from the steady-state; the time integrals in 
Eqs. (B1-B2) can then be split in two parts: 

Tr = A - t/£{C(r, t)} = J dt' J dr' 2 G(xr(f), xr 0 (0; * ~ 0 
+ £ dt' J dr' [2 -I- C(r\ *')] G(x r (t), x r '(f'); t - t') 

— 2A f dt' f dr' [ dz'G{xr(t),x';t - t') v 0 (x') ■V' J .To(x') 

J — 00 J o(r # ) 

-2A f dt' [ dr' f dz'G(xr(t), x'; t — t') v(x\t') •V z T(x\t'), (B3) 
Jto J >,f) 

V*T(x, t ) = J' 0 dt'J dr' 2 V x G(x, x r ' 0 (0; * - 0 


+ dt' j dr' [2 + C(r',0 V r G(x,xr — t') 

/ to r roo 

dt' / dr' / dz'V x G(x, x';t- t') v 0 (x') • V' x To(x') 

~oo J J<Z o(r') 

— 2A / dt' [ dr' [ dz'V x G(x,x'-,t-t')v{x',t')-V x T(x',t'). {BA) 

Jto j J C(r',l') 


Let t = iff — t' and 


xr(r,< n ) = x(r,C(r,t„)) = x„(r), n = 0,1,...,. 


The first integral in Eq. (B3) becomes 

Ii=f dr / <fr'2G(xN,x 0 ';r) 

J tp; — <o ** 

= j dr j dr' — exp < - 

Jtfi-to J A (ttt)* [ 


-|x N - x 0 ' 


- r - ( N + C' 0 > , 


where 


|xN-x n '| = |r-r'| + (Cn -Cn) > n = 0,1,..., JV, 


C(r',<«) = C 


Then 


|x N -x 0 '| j |x N -Xo'| 

lo = — — = — and ujn = . 

2x/T 2%/jVA? 


i x = r ^ exP - 

J 7r!|x N -xo'| Jo \ 


|x N - Xp' 
4a; 2 


f j r / ex P {~Qv( r ) + Co ( r> ) } 
J Air |x n - Xo'| 

x [exp |xn - X 0 '| erf (lo + - 


|x N - x 0 ' 


+ exp {- |x N - x 0 '|} erf ( lo - 


jxN - Xq > 

2 co 
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Finally 


h 


dr 


, exp{-(A’(r) + Co , (r / )} 


I m 4?r |xn — xo' 

* erf te^ +v7Fr7 ° 


— 2 sinh |xn — x 0 '| + exp |xn — x 0 ' 
+ exp {— |xn - xo'|} erf 


(Bo) 


>N -Xp , | — 1 1 

' / " °Jj 


The second integral in Eq. (B3) becomes: 


h = J drjdr' 2-f('(r G(x N ,x n ';T) 

~L dT I' 


rt N -t 0 /■ 2 + C(r',tft - t) 

dr / dr' j - 

'o J (47 tt) 2 


— l r — r'| 2 + [CA r ( r ) ~ CWitN ~ t) + 2r 

At 


We now discretize the time domain [0, t pj — to] in N intervals [(N — n) At, (N — n + 
1)A<], n = 1,2,..., JV. The value n = 0 corresponds to r = NAt = tpj, or, t' = t 0 
(initial time), whereas n=N corresponds to r = 0, or, t' = (final time.) Quan- 
tities, such as the interface location ((r, t), the temperature T(x, f), and the flow 
velocity v(x, <), are considered constant over each interval and equal to the values 
corresponding to the end of the time interval. The integral now becomes a sum of 
N integrals and the time integration can be performed analytically: 


N f(N—n+l)At t 

= V / dr I dr' 

= 5Z f dr ' 2 + £«(*■')] exp{-CA/(r) + Cn(r / )} x 

n= 1 ^ 


2 + Cn(r') 

V-exp 

(47Tt) 2 


r - r'l 2 + [Q/(r) - Cn(r') + 2 rf 


At 


k exp l _ ^ f (fa -a' _ r 

J(N- 


(N-n)At (47 Tt)’ 


4t 


Again, let 


|XN“X n '| j |X N -X n '| 

u ; = 1 7 = and , — =. 

2v^ 2 ^/(JV - n)Ai 


Then 


N r r 1 

h = £ / rfr ' [ 2 + Cn(r') 

n = l ^ 


exp{-CA/(r) + Cn(r')} 
27T f |x N - x n '| 


X 


r N " B+1 , / 2 |XN-X n '| 2 l 

l 4 ^ 2 J 


The integration in finally yields 

N 


l2 = tf *' [ 2 + <»(Ol 

n=l ^ 


nl ex P{-CN(r) + Cn(r')} 


8 tt |x n - x n ' 
x { exp |xn — x n '| I erf 


— erf 

+ exp {- |x N - Xn'l} ( erf 
— erf 


^N-V|. + V (w _„ )At 

2 y/(N - n)A t 

|xn-x„'| +V (JV _ n + 1)A< 
2 % /(iV-n+ 1)A< 

|xn ~ x "' 1 y ( jv-„jAt 

2-/(Af-n)A( 


(B6) 


; |XN - X "'' -V (AT-n + l)A t 

2y/{N — n + l)At 

We now follow the same procedure for the evaluation of the domain integrals 
in Eq. (B3): 

roc r roc 

= — 2A / dr dr' dz'G(xr(t),x';t — t') Vo(x') -X7' x To{x') 

JtN—ta J J Cn(r') 


f tN —to 

roo 


r oo /■> /»oo 

= — A I dr I dr' I 

Jttf — to J J Ca(r' 


dz' 


a ex P 

Co(r') 4(7 Tt) 2 


[x N - X 1 

4r 


»i 2 


-t-Cn + 2' Uo(x')-V;To(x'). 


Let 


xn-x'I j |xn — x'| 

u> = — 7 — — — and l>n = 


Then 


J 3 = -A / dr' r dz 
J dc 0 (r') 


2y/T 


,exp{-Cw(r) + z’} 


2y/NAt' 


7T2 |XN — X' 


-^vo(x') -V' To(x') 


P j (2 I x n-x'| 2 \ 

x jf ^exp|-u. j. 





The integration in u> gives 

~~ [-Mr) + *'} 

uz 

'Co(r') 


I 3 = -\[dr'f dz'^ 
J J Co(r') 4 * 

i- 


4tt |x n - x'| 
x ^ — 2 sinh |xn — x'l + exp {|xn — x'|} erf 


v 0 (x , )-v;to(x') 
|x n - x'l 


[2 x/tw — to 


+ vtjsr — to 


+ exp {— |xn —x'l} erf 


xn - x 


L2V<n - <0 

The second domain integral in Eq. (B3) becomes with r = — t': 


— y/tN ~ to }• (-B7) 


ftN—t 0 t tOO 

I 4 = — 2A / dr dr' dz'G(xr(<), x'; < — t') v(x', t') • V' I T(x', <') 

./o J J C(r' ,<') 

rttf—to r too 

— —A / dr dr' / 

« 7 o 4 ( 7 rr) 


^-yexp | - |XN i _. X ' 1 - - r - Cn + 2 ' } v(x'/) -W.C). 


We now apply the time discretization (see Eq. (B6)) and we re-introduce the 
variable transformation 


|xn — x'| |x N -x'| 

to = — r — — and u>N-n — 


2 y/r 


2 y/{N — n)A< 


Then 


A = -A £ / *' r dz ' exp ( ~ (w(r) + - ?- } v.(x') • V' T„(x') 

•'Cn(r') n 2 | X N — x'l 


XN - X'| 
x / du exp < — cj 2 — 

JuJN-n ( 


|x N ~ x 

4u> 2 


/ 1 2 


The integration in u> gives: 


, exp {— Gv(r) + 2 '} 


/ 4 = - a E / rfr ' r dz ' , SiVV * ; -v' x T n (x') 

An(r') 47r|x N -X'| 


x ^ exp {|xn — x'l) j erf 
— erf 


x N - x' 


2v/(A r - n + 1)A< 
+ exp {- |x N - x'|} ferf 


1 / + v(N — n)At 

[2-\/(iV — n)At 

+v / (W-n + l)A< 

A'-n + DAi 


(B8) 


|XN ~ X-1 v/ (JV - n)A< 

2>/(Ar-n)A< 


— erf 


x N - x r 


2^(A 7 — n + 1)A t 


- y/(N-n+l)bt 
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The discretized (in time) equivalent of Eq. (B3) is now 


T r = A - vfC{(( r, t)} = h + I 2 + I 3 + I 4 (B 9) 


where the integrals (I 4 — I 4 ) are shown in Eqs. (B5-B8). The calculation of the 
interface shape at time t = tpj, as represented by Cn(f) (or, xn,) requires the prior 
knowledge of the temperature gradient field in the interior of the domain, V x T(x, t). 
We shall use Eq. (B4) to evaluate the temperature gradient; our methodology will 
be similar to that applied to Eq. (B3). 

We again discretize the time domain in N intervals and the interface location 
C(r, t), the temperature T(x,t), and the flow velocity v(x, t) are considered con- 
stant over each interval. However, now these are considered equal to the values 
corresponding to the beginning of the time interval. This allows us to calculate 
the temperature gradient at time f.v (LHS) from quantities already calculated in 
previous time steps (RHS). •••■* 

The first integral in Eq. (B4) becomes, after the transformation r = t/v — 


Hi = f dr ( dr'2 V x G(x, Xc/;t) 

J t /v —to J 


= f dr f dr' — ^ 

Jifl—to j 7 T 1 


4(7Tt) 


rVx 

2 


exp 


-lx - X o 
4 r 


m2 


- T - Z + Co 


Instead of performing the differentiation in the spatial variables first and then inte- 
grating in time, we use the result of the time integration in I\ (see Eq. (B5)) and 
then differentiate (we first have to replace Xn with x since the (H) integrals refer 
to interior points): 

* - Co 


H i = 


M[ 


-A + (B - A) 


|x-Xo' 


l z + 


(B-A) 




(BIO) 



where 


^ - exp {-z + <„>')) j , , |erf 
477 |x — Xq'I [ 


+ exp{-|x-x 0 '|}erf J* = J - ~ *o 

zvz/v — to 


— 2sinh |x — x 0 '| j 


and 


B = -^exp{-z + <o( r ')} 


—2 cosh |x — X(/| 


4- exp |x — Xo'| erf 




xerf 




2y/tN ~~ h 


-exp{-|x-x 0 '|} 

M 2 


2 | ix-xo', 

+ / , . exp { -777 7-7-<N + <0 


\Zn(tN - t 0 ) 1 4(<a? - to) 


The second integral in Eq. (B4) is: 


H 2 = J'J °drj dr' [2 + C(r',0] V x G(x,x n ';r) 

-L dr J 


rt N -t 0 r 2 + C(r',<^-r) 

dr / dr' — 3 i 

(47Tt) 7 


x V, 


-|r - r'[ 2 + [z - C(r', t/y - r) + 2r] ii 
4r 


exp 


We now discretize in time and use Eq. (B6); again, we change Xn to x since x € ft. 
The crucial difference from I 2 is that here the quantities in the RHS are all known 
from previous steps: 


H 2 = ^ J dr' 2 + ( n _i(r')] ||-C„_i + 


(5 n -l — C„- 1) 


Z Cn — 1 


i , + (D n _i-C n -l) 


r — r cos <p 



ir r • (511) 


where 


exp {—z + 

'-'H — 1 o | /I X 

8tt |x - x n _i'| 

X {exp |x - *„-r'l (erf [ + VW^jKi 


|X- X n _i' 


— erf — . , = ===== + y (N — re + 1)A t 
2v/(A r -n + l)At VV ' 


+ exp {— |x — x n _i/|} (erf ~ VU* - n)At 


X - X n -l 


— erf — , == — W(N — n + 1)A t 

2y/(N-n+ TjKi VV ' 


-Dn -1 = exp {-2 + Cn-i(r')} X 

07 r 

X {exp |x - x„_r'| (erf [ A==^ + V(N - n)At 


|X- x n _i' 


-erf — .1 1 ■ + s /(N-n + l)At 

2 % /(iV-n + 1)A< V ' 


— exp {— |x — Xn-i'l} (erf ~ VW^t 


X- X n _i' 


+ - n)At eXP { ~ "S ' (N ~ n)A< } 

2 = exp ( Ix-Xn-i'l _ (jV _ n+1)A Al 

y/n{N — re + l)At 4(N — n + l)At jj 

The same procedure is followed for the domain integrals in Eq. (B4). The first 
domain integral becomes for r = tpj — t': 

H 3 = — 2A C dr f dr' f°° dz'V x G{x,x']t - t') v 0 (x') -V^x') 

JtN—tn J ‘'Co(r') 


1 


= -A / dr f dr' [ dz' — i-rvo(x') -V^x') 

Jtfj-to J J Co(r') 4(xr) 2 


x V x 


X - X , 

exp < t — z z 

At 


We now use the result of the time integration in I 3 (see Eq. (B7)), after we replace 
Xn with x, since x 6 fi: 


H 3 = -A / dr' /°° <fz'v 0 (x') • V z T 0 (x') 
7 J Co(r') 


-F + (F- F) 


(B12) 


z — z 


X — X' 


i* + (F-F) 


r — r 


cos</> \ . 


x — X' 



where 


JE = 


exp{- 


47T |x 


i+£l r 

- x'l l 


exp |x — x'| erf 


|x-x y l 

2\/F/v— 


+ — <0 


+ exp {— |x — x'j} erf 


x - x 


2 y/tw — to 


£ 


— 2 sinh |x 


- x,| l 


and 


F = 7 -exp {-2 + z'} 
47 T 


—2 cosh |x — x'| 


+ exp |x — x'| erf 


x — x 


2 jy 1 0 


+ y/tN ~ to 


xerf 


|x~ x'| 

L2>/<n ~ <0 


— \An — to 


- exp {— |x- x'l) 
-m2 


+ 


\A(<at - to) P | A(tN-t 0 ) 


X — X 


— <7V + ^0 


11 


The second domain integral in Eq. (B3) becomes for r = — t 


H 4 = -2A [ " ° dr f dr' f dz'V x G{x , x'; t - t') v(x', t') • V' t T(x\ i') 

7o 7 7c(r',t') 

f*N-to t f<X> J I 

= -A / dr dr' / ■ - v(x , ,< > ) -V^T(x',t') 

7o 7 J<;(r',t') 7T 2 |x - x' 


x V, 


exp ^ - ■ * ^ r - 2 + 2 ' 
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We now discretize in time and use Eq. (B8) (after we change xn to x, since x 6 fl). 
Once again, the basic difference from I 4 is that here the quantities in the RHS are 
all known from previous steps: 

N 


H 4 = [ dr ' [ dz 1 v n _ 1 (x') •V' I T n _i(x') 

n=l ^ •'Cn-i(r') 


— Qn - 1 + 


+ [Rn - 1 “ Qn- 1) 


Z — 2 


9 \ 


x — X' 


1* + {Rn- 1 “ Q n — 1 ) 


r — r cos <p \ . 

lx - xf 


(B 13) 


where 


_ exp {— z ■+ z'} .. 

t*Tn — 1 — .1 . 1 X 


47T X — x' 


x < exp |x — x'| ( erf 


X- X' 


[2V(A r -n)A< 

— erf 


+ V(7V - n)A< 


X — X 


2 V /(A T - n + l)At 


+ \/(iV - n + l)At 


+ exp {— |x — x'|} erf 


x — x 


[2 x /(A^-n)A< 
— erf 


- y/(N - n)At 


x — x 


2y/{N-n+ 1)A< 


— y/(N — n+ l)Af 


and 


{ 


Rn- 1 = — exp{-Z + 2 '} X 

47T 


x < exp {|x — x^} I erf 


x — x 


[2 x /(iV-n)A< 


+ y/(N - n)A t 


— exp {— |x — x'|} ( erf 


— erf 


x - x 


x-x' 


\_2y/(N — n)A< 


— erf 


[2>/(7V — n - h T)&i 
\ 

— y/{N — n)Ai 


+ - n -(- 1)A< 


x-x' 


[2y/{N - n + \)At 
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+ 


: exp 


lx - x'l 2 


y/n(N - n)At " r \ 4 (A 7 - n)At 

2 


- (N -n) At 


\/n{N - n + IjAt CXP \ 4(A 7 - n + l)At 


x — x 


— (N — n + l)Ai 


So the discretized equivalent of Eq. (B4) is now: 


V r T(x, t) = H 1 + H 2 + H 3 +H i 


(B14) 


with the integrals (Hi — H 4 ) shown in Eqs. (B9-B12.) 
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APPENDIX C: 


SPATIAL DISCRETIZATION 


This appendix contains the discretized equivalent of the spatial integrals de- 
rived in Appendix B. The coordinate transformations and the quadrature formulae 
that were employed can be found in Sec. 3.1 .b. 


N r A 7 *, 


'■-ZZ 


pTT 


2N N 

i=i fc=i ZiVri 'v’ 


tan 


2i — 1 7T 
2N~ r . 2 


ex P{— CN( r j) + Co(r,)} 


Y (%/<*«, A 7 + Wi + V’ifc.jv) (\/ Q 'i,N - A + ^ifc,7v) 


-2 sinh( w) + exp{ } erf 


. + \/aTa? 

[2\/TTKt 


+ exp{-V’,jfc,iv} erf 


[2y/7TKi 


■, (Ci) 


where 


Q'z,A 7 = r? + rj + [Cjv(tj) - Co(r,)] 2 and 0, - 2r,ry, 

2 A- — 1 7T 


V’iifc.TV = \A*i,N ~~0i+ \/ a i,N + Pi ~ y/ &i,N — fii 


cos 


2JV v 2 


A 7 N r AV 

*-zzz 


p7T 


i fi 4A7 '-^ 

n=l i=l *=1 ^ 


2 + Cn(r<) 


tan 


x 


/ 2i — 1 7T 

V 2N~r 2. 
ex p{— C n(^j) + Cwfo)} 


^/(\A*i,./V-n + ft + ^iJfc.AT-n) (\A*i,W-n “ ft + V’ifc.N-n) 


exp{V>;jt,N-n} 


erf ( — 

l 2y/(A — n) At vv ' 


X 



+ exp{^’iA,A'-n} 


- erf ( — ^ik,N-n _ -n+l)At 

\2y/(N-n + l)At VV ' y 

-- f ( 2 ^ , ::~; 1)A , -^-’ i+i)At )]}- (C2) 


where 


<*.,N-n = r? + r? + (Cn(^) - Cn(ri)] 2 and ft = 2 r^ry, 


*Pik,N-n = \Z<*i,N ~ Pi + \/ tti.JV-n + A ~ y^N^n~Wi 


COS 


2 A: — 1 7r 
2N V 2 


Av Av A r 

^3 = -aEEE 


07^ [ 2i — 1 7T \ /" 2m — 1 7T 

— tan I —m — — I tan 


jri ^ 4(1 + \)N r N v N z "" V 2 iV r 2 y-"V 2 JV, 2 

e xp{-CN(r>) + Co(ri) - 57^} 


x/CA im,N + Pi + V’i mA-./V MV ®im y N - ft + V’imfc.A?) 

-2sinh(V>i mfeiN ) + exp{j/w,A r }erf 


%££= + a/]va7 

2\/jVA7 


(C3) 


+ exp{-7/>i m A,N}erf 


. -\/NAi 

[2VNAi 


X 


v.(r' = r„ 2 ' = Co(r.)-^ly)j 


ViT 0 (/= rj , 2 ' = C.M-|%) 


where 


2 1 2 1 

Q'im,N = r t - + rj + 


Cw(rj) - Co(r.) + 


ln«m) 


1 2 


2(1 + A) J 1 


Pi = 2r,r 


J » 


1pimk,N — y/oiim,N Pi d" \Z®im,N ” Pi \f &im,N Pi 
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cos 


f 2k — 1 7T 

V 2JV v 2 



N N r A v- . N z 

*~*EEEE 




2i — 1 7r \ ( 2m — 1 7r 

tan — — — -- tan 


^ ^ ^ 4(1 + X)N r N v N z — V 2A r r 2 J "" V 2iV z 2 

^ ex P{ — CA'( r >) + Cn( r i) ~ 2(1 + A) } 

\[W ^‘m,A T - r? + A + ^im^,A T -n) ( \/ ^jm,A ? -n Pi "I" ^imk ,A T — rj) 


X < exp{V>, m fc,N-n} 


+ exp{V>, m *,N- n } 


f/ + >/(w-")aA 

-n)At V ' ' / 

- erf ( 7 =^= == + x/(iV-n + 1 )aA 

\2 % /(iV-n + l)Ai VV 7 

erf ( ■ ti = fc ’ N ~ - ^ = - >/(JV -n) Ai ) ( 

V2v/(JV-n)Ai 7 


- erf ( 7 Vw ’ N ~ n V(JV-n + l)A< 

\2^/(iV-n + l)At 


X 


(r' = ri - 2 ' = <"(r<)-^TIj) J 




where 


2 i 2 | 

a im,N — n = r, + ry + 


CA’( r >) ~ Cn(^) + 


HU) 

2(1 + A) J 


Pi = 2r,ry, 

1pimk,N — n = \f ® i m , N — n Pi 

y/Gir. ri^N — n~^Pi — n Pi 


COS 


2 /c — 1 7 T 

2A^ 2 


The vector integrals (i?i — #4), that correspond to the gradient of the temper- 
ature field in the interior of the domain, are dotted into the flow velocity 
to give the following discretized equivalents: 


N r N v 


*>-EE 


P 7 T 


, , 2N r N v 

H = 1 v—\ v 


tan 


( 2 M~ 1 *) 

V 2 N r 2 J 



ex P{~Ov( r «) + Co(r M ) + } 

\J (\/7m,A' + /?/, + # fiv,N ) (\/7/i,A' - Pn + 0 

in(^m) \ 


X < 


-A + 


5- .4 


B - A f (7 »,n - 0j„,w) 

w r ' 2 r 


+tv 1 


; (C5) 


where 


7/1, N = r 2 + rf + 


Cn(h) - Co(r„) - 


HU) 


9 


liv,N — 


2(1 + A) J ’ 

= 2 r„r,-, 

%/ 7/i, N - l^ft + [\/twn -TS”- %/tWv - ^ cos 2 | j 


The quantities A, A, and jB are given by 


A = < expj^^,^} erf 


0 


+ \/iV At 


L 2 v/]VA 7 


+ exp {~9 fll/t N ) erf 


0 


^ y/NAt 


[2^NAt 


- 2sinh(0 MI/>N ) 


A = 7— — I exp-^^A'} erf 
P/u/,N l 

+ exp{-O fil/iN ) erf 


v / 7VAt 

L 2 \/]va 7 


W ./^T 7 
2\fWKt 


— 2 sinh (0 




and 


B = 


2 cosh (^/ii/, a/ ) ~h €xp{9/xi/ t jv } erf 


9 


'^, + VNAt 


[ 2 \^NAt 


- exp {- 0 „„,Ar} 


x erf 


0 „ 


. p^=-VNAt 
[2VNAi 


H . exp ^ — N At 

V^C WAt [ l 4(iVAt) 


The velocity components v r and are evaluated at r = r,- and 2 = CA ? ( r «) — 2(1+™) 



N Nr AV 


"* = EEE 

a — 1 /j =1 i/=l 


pit 


4 N r N v 


tan 


2/1 — 1 7T 

~ 2 N 72 


2 + Ctr — 1 ( r #/ ) 


exp{— Cjv(r,) + C<t-i(*>) + 5Tifer} 


x < d, 


,A r — <7+1 ”1” “I" <7+1 ) <r+l ) 


+ ,v e 


Do - 1 — C(T-1 


/U7,N-<r+l 


rv - 


(7 /*,n-«m-i - 0j 


liv,N — < 7 + 1 ; 


2ri 


(C6) 


where 


7 ft, A?— 17+1 = rj + r? + 


C N{ r i) C<7 — l( r ft) 


M€m) 
2(1 + A) 


fin — 2r M rj, 

^ftt/,A r — tr+l = \/ 7 /i,A?-<t+1 — fin + 

y/^fn^N— tr+i 4" ^/t \/ 7/t, A^— <7+1 ~ /^/t 

The quantities C c -i, Co- i, and are given by 


cos 


2z/ — 1 7T 

2iV v 2 


CV-i = < exp{^ t|/i A/_ <r+ j}( erf 


0 


r A, ~ g+ 1 +\/(iv-a)A< 

2 ^(JV - a)A< 


+ exp{-^„ )W _ (T+1 } erf 


— erf 
0 


- fefatV — + ^/(jV - a + 1)A t 

2 y/(N - a + l)At 1 


4£=£= — - V(7V - g)A< 
2 y/(N-r*)At 


— erf 


. ^(AT-ff + lJAt 

[ 2 V(JV - ff + l)At 


— f 

t/,N-er+l ^ 


Cff-i = 7 1 \ exp{0 M( ,,A/_ ff+1 } erf 

t'/ii/.N-i 


6 


w-»+i ^ (iV _ 

2 ^(JV ~ <r)Ai ’ 


2 % /(7 V-< T + i)At VV 1 - 


— erf 



+ exp {-0#i«/,A'-«h-i} erf 


7===-VW^Wt 


_e_ 

2 y/{N-a)At 

e 


— erf 


. r£z°±' -v/QV-a + ^Af 

2 % /(A r -a+ 1)A< 


and 


^<7 — 1 — \ €Xp A 7 — <7+1 } ( erf 


e, 


. f v ' N -* +1 + y (A r _ ct)A ; 

[2^/(A r - a)A* 


exp { <7-fi } I erf 


— erf 
6 , 


a 

2VP 


b ^g± L = + V(N - <r + l)At j 

\ r — <r + i)a< Jy 


g+1 — _ Jjjj _ g)A< 
2 v/(A^-<7)A < 


— erf 


fefe±i== - V(F37 TT)a7 


0 2 . 


\Ar(A r - <r)At CXP \ 4(AT — a)A< 

2 

: exp 


ftt'. A 7 ~ <7 4-1 


2 y/(N — a + 1)A< 
- (A' - a)A< 


Z}2 


y'V(A r — <7 + 1)A£ ( 4(A r — <7 + l)A£ 


— (Af — a + l)Af > >. 


The velocity components v T and v z are evaluated at r = and z = C/v( r j) — 20+A 7 


£t __\V'V'V' E-1 t / 2p - 1 7T \ / 2/ - 1 7T 

3 ~ U 4 (! + A)JWMr. v 2Ar r 2 J V 2iV z 2 


exp{-CAr(r,) + Co(r M ) + 


x < Vz 




-E+ A-5(< K ( ri )-<o(r„)- 


Ji'Jlpt.N +011 + Bplf,! v) ( — 


+ Ur 


2(1 -t- A) 
F — E 

0tilv,N 


) 


(t/W.A' - ^/„,A') 


2r, 


>, (C7) 
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where 


V 


TftiLN = r 2 + r? + 


CaK**.) - CoM - 


l°(U/fr) 

2(1 + A) 


1 2 


/?/i = 2 r ll r i , 

Olitv,N = \/ 7 #*/, AT -/?/!+ \/ + Wn - V 7 /i/,N - Pp 
The quantities F, F, and F are given by 


cos 


2 */- 

2iV, 


- 1 7r\ 


= jexp^j,,,^} 


erf 


0 


i^= + VNAt 


[2VNAt 

+ exp erf J;^==== - ViV At - 2sinh(^ M ^ jN )|, 

£ = ^ “ f [s^B + 

- 2sinh(^/ V| jv)|, 


+ exp {—Out,,'!*} erf ^7== - s/TTKi 


and 


F = 


—2 cosh^/^w) + exp-ffl^ArJerf 


[2VNAt 


x erf 




- VN At 


+ VNAt 


el, 


— exp {— 0 M fj/,w} 


, 2 I tj A , 

+ ^fBj eXP \-4(lVAi)- ArAt 


Un/FaT 

The velocity components u r and t>- are evaluated at r = r* and z = £jv( r ») — 2 (i+a) • 


i. 


N N r A% w, 

**~*E£EE 


p7r 


tan(H^-Vn(— l) 

\ 2 N r 2 J \ 2 N Z 2 J 


XZZiiZ* 1 + W^JV. 

. exp{ -Cn (r,Q + +ffl } 

^(\/7j</,A/-<r+l +^7+^/i/i/,N-<r+l) (\/7 /i»,N-<t+ 1 ~ /?#x + ^Jt/,N-<T+1 ) 


X < V, 


— Qct-1 + 


Re — 1 Q<t-1 /a / \ > / \ l n (£m/£/)\ 

^ r (a.(r.)-c.- 1 (r„)- 1(TTTr j 


4- V 


e^iv 

Qa-1 ~ Re-1 


r i - 


(lnt,N-e+l ~ e‘j llViN _ a+1 'j 
2ri 


(C 8) 
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where 


7/i/,A ! -er+l — + >7 + 


Ca r ( r t) Ca— i( r /i) 


MUM 

2(1 + A) J 


T 2 


P* = 2r„ri, 


6 


<r+l — 


\/ 7 M *,A r -(7+l - /?*/ + [\/7/W,N-<7+l T^T- >/ 7 m/,N-<t+1 - S] COS 2 — 2”^ 


The quantities Ccr-i, and -R^-i are given by 




Qtr — 1 — 


|exp{« p , (erf + V(* - ^ 


— erf 


yu,N-*+ } + y /(N-a + l)At 

2-v/(JV -cr + i)A< 


+ exp {-^/^N-a+i} ( erf 


BrW-e+l 
2^/(iV -<r)A< 


— erf 


<W-. + l -._^ (JV _ g+1 ) A< 

2\/{N — <7 -{- l)Atf 


0<7 — 1 — 


— cr+1 


ex p{^/i/i/,N-(T+l } erf 


^<i/,A r — g-fl 

2V(A r - a)A* 


+ >/(JV - a)At 


— erf 


W-.+1 ^ (A f_ g+ l) A< 1) 

2^J{N — a + l)Ai V ' ’ \) 


+ exp { ( 7 +i > erf 


^,N- g +i ; _^tT7)a7 

2^/(^ -<7)A< 


— erf 


0 


2V(W - ff + l)Ai 


and 


I 


#< 7-1 = S exp{0 #l/l/iN _ (y+ i} erf 


- q)Ai 

2\/(N - <r)Ai! V ' 
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— erf 


2\/(A r -<T+ 1)A( V ' \) 


exp { 8 u . j_i } erf 


e ^,N-a+l = _^ (A r_ a)A< 
2y/{N - a)At 

6 


— erf 


^ti/.A 7 — <t+i _^ (A r_ a + 1)A< 

2y/{N — <7 + l)At 


, exp / — T7Tr — — TT (-W — (T')A^ 

•^/V (JV.— cr)At 4(iV — cr)At 


0 2 . 


y/n(N — (7 + l)At ^ 4(JV — <7 -f l)At 

The velocity components u r and u z are evaluated at r = r, and z = CAf( r «) ~ 2(1 +a) • 
The expressions Ci through Cs can now be combined to represent the dis- 
cretized equivalent of the integral evolution equations. 


1-1 


-(N-a + l)At 
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APPENDIX D 


Boundary integral formulation for the fluid velocity 

In this appendix we present the boundary integral formulation used to evaluate 
the velocity field v at each time step. The flow past the unperturbed paraboloid is 
represented in the entire fluid domain by the Oseen approximation to the Navier- 
Stokes equation. When the crystal-melt interface is perturbed from its paraboloidal 
shape, one has to correct the flow field so that the no-slip condition on the perturbed 
surface is satisfied. The flow velocity v can then be expressed as a sum of two 
contributions: 

a) the Oseen velocity, Vo, for the unperturbed paraboloid , which is known in ana- 
lytical form but does not satisfy the no-slip condition on the perturbed surface, 
and 

b) a “correction” that vanishes at infinity but is equal to the negative of the Oseen 
velocity on the surface. 

The total flow velocity then becomes 

v = v 0 + u (-D1) 

with 

u = — v 0 , on the interface (D2) 

and 

u — * 0, as x — > oo. (-03) 

The velocity u is only a local correction to the Oseen flow field in the vicinity of 
the surface; the asymptotic behavior of the total flow field at large distances from 


the solid is determined solely by the Oseen solution. Thus, for the case of very low 
Reynolds numbers present in dendritic growth, the flow field u is inertialess and 
can be described by the Stokes equation 


-VP + V 2 u = 0. 


m 


(As was discussed in Chapter 3, the temporal term of the transient Stokes equation 
can be neglected because, at very low Reynolds numbers, the velocity field responds 
instantaneously to changes in the interface shape). 

We solve the Stokes equation using the boundary integral formulation first 
applied by Youngren and Acrivos [1]. The differential equation is transformed into 
a boundary integral equation that relates the local surface stresses f to the velocity 


:( x)u(x)=±/ 


3 / (x-y)[(x-y)(x-y):n(y)u(y)] 


+ 7T” 


1 r [f(y) . (X - y) [(x - y) • f(y)] 


r L a xy 


dr, y € T, (D5) 


where n is the unit normal on the interface pointing toward the fluid domain £2, 


d xy = |x - y|, and 


_ / j, if x € T; 

\ 1, if x € £2 — 


The calculation of the interior velocity takes place in two steps. First, the surface 
stresses, f are calculated from Eq. (D5) with x € T (the surface velocities u(xr) 
are known from the boundary conditions, i.e. the Oseen flow velocity). Then, Eq. 
(D5) is used to evaluate the interior velocities for selected points x 6 £2 — T. 
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For axisymmetric flows, the vector equation Eq. (D5) reduces to two scalar 
equations. In a polar cylindrical coordinate system, where x = (r x ,ip x ,z x ) and 


y = ( r y > ^Pyi z y)> Eq. (D5) can be decomposed to 






z x C y ( r x cos $ r y) 


x [u r (y)(r* cos 0 - r y ) + (** - Cy)«*(y)l 


= hl de l dr ' r »v 1 + 


OF 


(y) cos 0 i r x —r y cos 0 

+ W, 


x [/r(y)(r*costf-r y ) + (« x -C»)/*(y)]>, (D6) 


)■ 


and 


:(x) “' (x) -^r^r% 7 ( ^-« 


0 4y 


^x Cy ( r x cos ^ r y) 


x [u r (y)(r* cos 0 - r y ) + ( 2 * - C»)«*(y)] 


“8 ll i9 l dr » r »^ I+ (*') {4^ 


/z(y) . Cy 
d 3 

“xy 


X [/r(y)(r.costf-r y ) + (z*-C(y))/*(y)] (-D7) 


where z y = £y = C( r y) represents the interface T, 0 = — <^j,, and 


dxj, = y/rl+^l- 2r * r y cos 0 + (** - Cy) 2 • 


We now discretize the interface T in N r ring- type elements. Each element Tj is 
associated with a “node” at r = rj. Over each such element the surface velocities 
and the surface stresses are considered constant and are represented by their values 
at the corresponding node. The surface integrals in Eqs. (D6-D7) can now be 
written as the sums of the integrals over each surface element Tj. For x E T, Eqs. 


(D6-D7) become 


u r( r >) - 7T~ [ dQ f - r cos 9) [(,• - £(r) - C'( r )( r »' cos 0 - r)] 

Jo Jo a zy 


x [u r (r)(r,- cos 9 — r ) + (C, - C{r))u z (r)] 




xy 

2n 


cos 9 + 


r-j — r cos 9 


d \ v 


(D 8 ) 


+ E M r i) l M l ^-\A+(C'M) J M - r cos 9) (C. - CM) 


and 


u *( r ^-'ixJ 0 dd J 0 ^ j "^ <_ ^ r ))^‘“^ r ) _ ^ r ^ r<cos ^ _r ^ 

x [u r (r)(r, cos 9 - r) + (C,- - C( r )) u *( r )] 
= XJ/>-( r j)^ d0 J r Jplf/ 1 + (C'(r )) 2 ( r < cos * - r) (Cl - C(r)) (-D 9 ) 


dr r 


+ 


E/.-M) f <* / 4^/TT«W ? [i + 

J—J »/0 ^ 7ra xy 


(g z^ir 
4, 


where ri = r x ,r y = r, z x = C(r.) = £'(r) = and 


d XJ , = yi r? + r 2 - 2r^r cos 9 + (£,• - C( 7 ’)) 2 - 


To calculate the left-hand-side of Eqs. (D8-D9) one needs to obtain the boundary 
condition for the velocity u. From Eq. (D2) 


u(x r ) = -v 0 (x r ). 


But v 0 is by definition equal to the Oseen velocity that satisfies Eq. (2.13) and 
is presented in analytical form in Eq. (2.14). This velocity is non-zero on the 
perturbed interface z = £(r), since it satisfies the no-slip condition on the surface 



of the unperturbed paraboloid z = Co( r )- Equations (2.14) and (D2) are thus used 
to calculate the boundary values of u as explicit functions of the interface shape 
z — ((r). (The Oseen velocity in Eq. (2.14) is expressed in paraboloidal coordinates 
and a simple transformation to the cylindrical coordinate system is required). The 
expressions in the left-hand-side of Eqs. (D8-D9) are then calculated using a simple 
Gauss-Legendre quadrature. 

The components of the local surface stress vector, f r and f z , are considered 
constant over each element Ty and are represented by the nodal vectors 


tr = [/r(ri),/ r (r 2 ),...,/r(rN r )], 
t z = [/r(ri),/,(r 2 ),...,/,(rjv p )]. 


The choice of the nodal points (r t ) in this problem is dictated by the discretization 
of the boundary integrals in the integral evolution equation for the interface shape. 
Equations (D8-D9) form a (2A r r x 2N r ) system of linear algebraic equations for the 
local surface stresses (written in compact notation) 


Ai A 2 \ f t r 

A3 A4 J \ t x 


(D 10 ) 


where the matrix element [A^ij], k = 1,2, 3, 4, represents the integral over the 
element Tj with r x = r,-. 

The integrals Ak,ij can be easily calculated for i ^ j using a simple Gauss- 
Legendre quadrature. For i — j, the integrals Atjj are singular because the inte- 
gration interval includes the field point r,- and thus d xy in the denominator becomes 
zero as r — ► Tj. Nevertheless, the improper integrals Ak,jj do exist as long as the 


surface T satisfies the Lyapunov condition at each node r = rj. It can be easily 
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shown that the integrands in each of the singular integrals behave as 1 /d xy in the 
vicinity of the singular point (r = rj). This weak singularity is integrable and the 
integration is carried out using the coordinate transformation 


where 


and 


The integral 


V = \Mr, r<) - 7 (r,r,)cos 0 ), 


a(r,rj) = r 2 + r 2 + (ty - C0*)V 


7 (r,r,-) = 2 r,r. 


= f d6 I dr 

Jo JTj d 


h(r,6) 


where h is a non-singular function, becomes 


/= I dr, / , iT , 4 ' l(r ’’ ,) 

Jrj \/(\/ a ' + 7 - *?) (*? - >/« - 7 ) \/(\A* + 7 + »?) (*? + \Z<* - 7) 

We can now apply following Gaussian quadrature formula (based on the orthogonal 
Chebyshev polynomials of the first kind) [2] 


L 


My) 


Vr, 


a V / (y-°)( t -y) m =l 

(6-fa) , (6- 


M 

£ jvjMym) 


+ 


a) / 2m — 1 7r\ 
— cos — — — — 

\ M 2 ) 


to the integral I to get 


4tt y-v h(r, q m ) 

M “"i \A\A* + T+ ym)(ym + vte-7)’ 


(Till) 
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where 


\/a + 7 + y/a- - 7 \/a + 7 — \/a -.7 / 2 m — 1 7 r \ 

’’” = — + 2 C ° S \ ) 

We now calculate all the integrals Ak,ij and then invert the matrix Eq. (DIO) to 
obtain the local surface stresses at each node. 

In order to calculate the velocity at a point x = (r,-, 2 /), in the interior of the 
domain, we use Eqs. (D6-D7) for x € D — T: 


u r (ri,zi) = 

7-/ d 6 f ~ r cos 9) [2/ — C( r ) — C ( r )( r « cos ^ “ r )] 

Jq J o a~ y 

x [u r (r)(r,- cos 9 - r) + (2/ - C(r)) u*(r)j 

\ 7 . 

r ,• — r cos 9 


+ Y^fr(rj) f del + (C'(r )) 2 cos 9 +^—^- 

*/0 OTTrf xl/ L a xi 

+ iMn)f 0 dd J r g^fV 1 + (C'( r )) 2 ( r » — r cos 9) (zi — C( r )) (£ 12 ) 


and 


4tt 


r2ir roc 

d6 

Jo Jo 


dr r 

17 


u z (ri,zi) = 

( 2 / - C(r)) [zi - C(r) - C( r )( r i cos 9 — r)] 


x [tx r (r)(r,- cos 9 - r) + (2/ - C(r)) u*(r)j 
+ Y^fr( r i) J Q d9 ~jry/ 1 + (C'(r)) 2 (n cos 9 - r) (2, - C(r)) 


1 + 


jzi - C( r )) 2 


<*s y 


(£13) 


where 


d xy = \f rj + r 2 - 2r,r cos ^ + (2, - C(r)) 2 . 



Equations (D 12-13) do not contain any singular integrals (x 6 fi - T whereas 
y £ T) and the Stokes -correction velocities can be easily calculated using a Gauss- 
Legendre quadrature: The total velocity v, which appears in the interface evolution 
equations in Chapters 2 and 3, is then calculated by adding the unperturbed Oseen 
velocity v 0 , given by Eq. (2.14), to the Stokes correction u for the perturbed inter- 
face. The new velocity v satisfies the no-slip condition on the perturbed interface 
and reduces to the unperturbed Oseen velocity away from the interface. 
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APPENDIX E: 

Derivation of the Gibbs-Thompson relation 

This appendix describes the derivation of the Gibbs-Thompson thermodynamic 
equation, which determines the effect of capillarity on the temperature of a crystal- 
melt interface. The derivation presented here is similar to that of Delves [1]. 

Consider a very small section of the interface, over which the curvature is 
constant, and with small adjacent regions of the solid and liquid phases which 
are uniform in temperature and composition. The total extensive thermodynamic 
variables ( E , S,V, F = Helmholtz free energy) can be divided between solid, liquid, 
and interface. Assuming that the volume and number of particles allocated to the 
interface are both zero, the surface free energy is defined by 

F(total) = Fs + Fl + 7 • Area, (FI) 

where 7 is the solid-liquid interfacial tension. The equilibrium conditions at the in- 
terface are obtained by minimizing the free energy with respect to changes of energy 
and of particles between the two phases. This leads to the equality of temperature 
and chemical potentials of each component (usi, Hu, i labels component), which 
are defined per atom of component. A further minimization with respect to an 
infinitesimal change of volume of one phase, the total volume being constant, leads 
to the condition 

Ps-PL-y—z r s = t (- + -)- (B2) 
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Ps and Pi are the pressures at the curved interface; rq and r 2 are the principal 
radii of curvature and are positive if the interface is concave towards the solid. If 
the pressures Ps and Pi are nearly the same as the pressures at a flat interface, 


A P s = Ps~ Ps(flat); |AP S | < P S , (EZ) 


one may use an approximate formula (approximate because A Ps is a small but finite 
change) from the Gibbs-Duhem relations, which are thermodynamic identities for 
each phase, 


£ A usi N s , = V s A P s - S s AT, 
A fin Nu = V L A Pi - S L AT. 
At a flat interface Ps(flat) = Pi(flat), so from Eq. (E2) 


(EA) 


AP s -APl = 7 ( — + — ), (£5) 

Vn r 2 J 

and A /is, = Ann because the phases are always in equilibrium, and A fiSi is the 
change in chemical potential if the interface changes from flat to curved. 

If the system has one component only (i.e. solidification from a pure melt), 
only one variable, A Pi for example, needs to be specified to find the other three 
variables AT, Am and A Ps- Assuming that the liquid is under a constant pressure 
and that variations in Pi at the interface due to hydrostatic effects and convection 
are negligible compared to A Pi, then it is reasonable to set A Pi = 0 and all the 
pressure difference is taken up by the solid phase. The elimination of A/J. from Eqs. 


Tm[Sl/Nl — Ss/Ns) is the latent heat of solidification per atom, so if L is defined 
as the latent heat per unit volume of solid and AT is the temperature at a curved 
interface minus the value at a flat interface, 

f r = Tm + a T = T„ (l - Ik) , (El) 

where Tr is the temperature at the interface T and K. — (1 /rj -f l/r 2 ) is the local 
interface curvature. By defining T = (T — T<x,)cp/L, A = (Tm — T^Cp/L, and K, 
= K./ p, Eq. (E7) becomes 

T r = A - (d 0 / P )K, (£8) 

where do = 7 Tmc p /L 2 is a capillary length and p is the tip radius of the crystal 
dendrite. Equation (E8) represents the Gibbs-Thompson relation in the case of a 
pure melt and is used in Chapter 3 as the interface boundary condition for the 
temperature field. 
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Fig. 1. Measured tip velocity vs supercooling for succinonitrile. Note the deviation 
from the pure diffusion steady-state theories at low supercoolings (reprinted 
from Glicksman and Huang, Acta Metall. 29 (1981) 701). 

Fig. 2. Measured tip radius t>,s supercooling for succinonitrile. Note the deviation from 
the pure diffusion steady-state theories at low supercoolings (reprinted from 
Glicksman and Huang, Acta Metall. 29 (19S1) 701). 

Fig. 3. Effect of convection on the morphology of the interface. The distortion of 
the temperature field causes partial elimination of the sidebranches (reprinted 
from Glicksman and Huang, Proc. 3rd European Symp. on Material Sciences 
in Space, Grenoble, April 1979). 

Fig. 4. Ratio of measured Peclet number to Peclet number predicted by the Ivantsov 
theory. Convective effects cause a sharp transition to occur at A = 0.05 
(reprinted from Glicksman and Huang, Proc. 3rd European Symp. on Ma- 
terial Sciences in Space, Grenoble, April 1979). 

Fig. 5. Diagram of the coordinate system showing the orientation of the paraboloidal 
crystal and the flow. 

Fig. 6. The Peclet number as a function of supercooling for selected values of the ve- 
locity ratio. For these calculations the Prandtl number is 23.2 (succinonitrile). 
The Ivantsov solution corresponds to A = 0. 

Fig. 7. The perturbation tip velocity £(0, t) is plotted us time for different numbers N r 
of discretization points in the radial direction (u = 0.001, A = 0.0, Pe = 1.0, 
N v = 7, A t = 0.002). 
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Fig. 8. The perturbation tip velocity £(0,f) is plotted vs time for different numbers 
A 7 ^ of discretization points in the ^-direction (u = 0.001, A = 0.0, Pe = 1.0, 
N r = 12, At = 0.002). 

Fig. 9. The perturbation tip velocity C(0,i) is plotted vs time for different time steps 
At (i / = 0.001, A = 0.0, Pe = 1.0, N r = 12, N v = 7). 

Fig. 10. Interface destabilization via tip splitting. The interface shape is plotted vs the 
radial coordinate r for different times (u = A = 0, Pe = 1, A r r = 12, N v = 7, 
A t = 0.002). 

Fig. 11. The perturbation growth velocity £(r, <) I s plotted vs the radial coordinate for 
different times and for A = 0.0, v = 0.020, and Pe = 1. 

Fig. 12. The perturbation growth velocity £(r, <) is plotted vs the radial coordinate for 
different times and for A = 0.0, v = 0.001, and Pe = 1. 

Fig. 13. The perturbation growth velocity C( r ^) I s plotted vs the radial 'oordinate for 
different times and for A = 0.1, v = 0.001, and Pe = 1. 

Fig. 14. The perturbation growth velocity C( r i*) I s plotted vs the radial coordinate for 
different times and for A = 1.0, v = 0.001, and Pe = 1. 

Fig. 15. The perturbation velocity at the tip, £(0,i), is plotted vs time for different 
values of the surface tension parameter v (A = 0.0, Pe = 1.0). The inset, which 
is a blow-up of the v = 0.001 curve, shows that anisotropy has a negligible effect 
on the growth velocity. 

Fig, 16. The perturbation velocity at the tip, £(0, t), is plotted vs time for different flow 
strengths A. Its absolute value increases with A, causing a decrease in the total 
growth velocity, [2 -f £] (u = 0.001, Pe = 1.0). 
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Fig 1. Measured tip velocity vs supercooling for succinonitrile. Note the deviation 
from the pure diffusion steady-state theories at low supercoolings. (Reprinted 
from Glicksman and Huang, Acta Metall. 29 (1981) 701). 
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Fig. 2. Measured tip radius vs supercooling for succinonitrile. Note the deviation from 
the pure diffusion steady-state theories at low supercoolings. (Reprinted from 
Glicksman and Huang, Acta Met all. 29 (1981) 701). 



Fig 3. Effect of convection on the morphology of the interface. The distortion of the 
temperature field causes partial elimination of the sidebranches. (Reprinted 
from Glicksman and Huang, Proc. 3rd European Symp. on Material Sciences 
in Space, Grenoble, April 1979). 
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Fig. Ratio of measured Peclet number to Peclet number predicted by the Ivantsov 
theory. Convective edicts cause a sharp transition to occur at A = 0.05. 
(Reprinted from Glicksman and Huang, Proc. 3rd European Symp. on Mate- 
rial Sciences in Space, Grenoble, April 1979). 




Fig. 5. Diagram of the coordinate system 
crystal and the flow. 






Fig. 7. The perturbation tip velocity £(0,<) plotted vs time for different numbers TV, 
of discretization points in the radial direction (u — 0.001, A = 0.0, Pe — 1.0 
Ny = 7, At = 0.002). 



Fig. 8. The perturbation tip velocity £(0,<) is plotted vs time for different numbers 
Ny> of discretization points in the (^-direction ( v = 0.001, A = 0.0, Pe = 1.0, 
Nr = 12, A t = 0.002). 
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Fig. 10. Interface destabilization via tip splitting. The interface 
6hape is plotted v$ the radial coordinate r for different 
times (u = X = 0, Pe = 1, N r — 12, N v = 7, At = 
0 . 002 ). 
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Fig. 15. The perturbation velocity at the tip, <(0,t), is plotted vs time for different values of the surface tension 
parameter v (A = 0.0, Pe = 1.0). The inset, which is a blow-up of the v - 0.001 curve, shows that 
anisotropy has a negligible effect on the growth velocity, (sec Eq. (3.3) for the definition of c). 
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Pe = 1.0). Its absolute value increases with A, causing a decrease in the total growth velocity, [2 + C]. 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


C APPENDIX G: LISTING OF THE FORTRAN CODE C 
C C 
C THIS PROGRAM CALCULATES THE EVOLUTION OF THE DENDRITIC C 
C INTERFACE IN TIME AND SPACE UNDER THE INFLUENCE OF C 
C FLUID FLOW AND SURFACE TENSION. C 
C C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
IMPLICIT REAL* 8 (A-H.O-Z) 

REAL*8 WK( 1700) ,UGT(30 ,30,1000) ,UGTO (30 , 30) , Z (1000 , 30) , 

. ZD ( 1000 , 30) ,X(30) , ZO (30) ,R(30) ,CSI (30) ,CS2(30) ,FI (30) , ZDINIT(30) 

, ,ORD(30) ,C(30,3) ,VECRH0(3) ,ZETA(3) , 

. BPAR (4 ) , STRES ( 60 ) , FLOW (60,60), QUAD ( 30 ) , WKAREA ( 60 ) , LAM , NU 
COMMON / A/UGT , UGTO , Z , ZD , ZDINIT , CS I , CS2 , FI , ZO , LAM , PI , 

. DT , NU , AN I S , NF , NZ , N , I NDEX , IVAN 
. /B/C ,R ,ORD,BPAR 
./D/ SMALL, P 

. /C/QUAD, RJ,ZJ,R2, MM, NR, UPFLOW 
EQUIVALENCE (WK.FLOW) 

INTEGER NPL(IOOO) ,NPERT( 1000) 

LOGICAL IVAN, UPFLOW 

EXTERNAL DCADRE , FCN , FA , FB , FC , FD , RIGHT1 , RIGHT2 , EXPIN 
CALL XUFLOW(O) 

WRITE ( 6 , * ) ' NR , NF , NZ . MM , NTOTAL , DT , PECLET , LAMBDA , NU ' 

READ (5,*) NR, NF.NZ, MM, NTOTAL, DT,P, LAM, NU 
READ (15,*) UPFLOW, PR, ANIS, (ZDINIT(I) , 1=1, NR) 

WRITE (6,*) 'THE TYPE OF FLOW (UPFLOW?), THE PRANDTL NUMBER, 

.THE DIM/LESS SURFACE TENSION, AND THE ANISOTROPY ARE:' 

WRITE (6 ) * #V • * v * V */r Vr “.v Vr -Jr -,V ~V -Jr -Jr ‘V “V ~ V V? it i% i% *V * V Vc Vr *V Vr “V *V "Jr “V “V it “V "V *V iV “V ^ Vr Vr * 

WRITE (6 , * ) UPFLOW , PR , NU , AN I S 

WR I TE (6 *V ^ * ‘Jr * Wr "V -Jr *,V -Jr * V -'r -Jr -Jr -Jr -Jr 7 V ~ V Vr -Jr -Jr -Jr ‘Jr * V ‘Jr -Jr -Jr -Jr -Jr "Jr -Jr V* "Jr -'r “ V * V * V <V "Jr Vr * V “ V » y /V * V -Jr V * * 

WRITE(6 ,*) 'IF YOU WANT TO CHANGE SOMETHING, HALT EXECUTION 
. AND GO TO FILE 15 AS DEFINED IN YOUR EXEC FILE 
. (SOME DATA USEFUL TO THE COMPUTATION ARE STORED IN FILE 14)' 

READ (14,*) SMALL, ZEROPL,RMAX, NS IG,ITMAX, 

. (BPAR(N) ,N=1 ,4) , (NPL(N) ,N=1, NTOTAL) 

READ (17,*) PRT, (NPERT(N) ,N=1, NTOTAL) 

IF (LAM.EQ.O.) I VAN= . TRUE . 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C CALCULATION OF DELTA (IVANTSOV) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DELI VA=P*DEXP ( P ) *EXP I N ( P ) 

PI=3. 14159265358979323 
RE=P* LAM/PR 
R2=RE/2 . 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C GENERATION OF ARRAYS FOR VELOCITY, QUADRATURES, ETC 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DO 810 M= 1 , MM 

810 QUAD ( M ) =DCOS ( DFLOAT ( 2*M - 1 ) *P I / 2 . /DFLOAT ( MM ) ) 

DO 5 1=1, NR 
R ( I )=DSQRT ( - 2 . *P*DLOG 

. (DCOS (DFLOAT ( 2* I - 1 ) /DFLOAT ( 2*NR )*PI / 2 . ) ) ) 

Z0(I)=P/2.*(1 . -R(I )**2/P**2) 


DO 5 M=1 , NZ 

CSI (M)=DC0S(DFL0AT(2*M- l)/DFLOAT(2*NZ)*PI /2 . )**2 
CS2 (M)=DTAN (DFLOAT (2*M- l)/DFLOAT (2*NZ)*PI/2 . ) 

IF (IVAN) GO TO 5 

ZIM=ZO ( I ) -DLOG (CSI (M) )/2 . / ( 1 . +LAM) 

W= ( ZI M+DSQRT ( Z I M**2+R ( I )**2 ) ) / P 
S= ( - Z I M+DSQRT ( Z I M**2+R ( I ) **2 ) ) / P 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C FLOW DOWN 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

IF (UPFLOW) GO TO 51 
UGTO(I ,M)=-2 . / (S+W)* 

. ( (DEXP ( -R2 ) -DEXP ( -R2*W) ) /R2/EXPIN (R2 )+ 

. W* (EXPIN (R2*W) /EXPIN (R2 ) - 1 . ) ) 

. *W** ( - 1 . +PR*DEXP ( -R2) /EXPIN (R2 ) ) 

. *DEXP ( P* ( 1 . +LAM ) * ( 1 . - W )■ + 

. PR/EXPIN (R2)*( -EXPIN (R2)+EXPIN2(R2)+ 

. EXPIN (R2*W ) -EXPIN2 (R2*W ) ) ) 

GO TO 5 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C FLOW UP 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

51 UGTO (I,M)=2./( S+W )* 

. ( (DEXP ( -R2 ) -DEXP( -R2*W) )/R2/EXPIN(R2)+ 

. W* ( EXP I N ( R 2*W ) / EXP I N ( R2 ) - 1 . ) ) 

. *W** ( - 1 . -PR*DEXP ( -R2 ) /EXPIN (R2 ) ) 

. -''DEXP ( P* ( 1 . -LAM ) * ( 1 . -W) - 
. PR/EXPIN (R2)* ( -EXPIN (R2 )+EXPIN2 (R2 )+ 

. EXPIN (R2*W) -EXPIN2 (R2*W) ) ) 

5 CONTINUE 

R(NR+l)=R(NR)+(RMAX-R(NR))/8. 

ORD(NR+l)= 

. P / 2 . * ( 1 - R ( NR+ 1 ) ** 2 / P** 2 ) 

R(NR+2)=R(NR)+(RMAX-R(NR))/3. 

ORD(NR+2)= 


. P/2 . * ( 1 -R (NR+2 )**2/P**2 ) 

R(NR+3)=R(NR)+(RMAX-R(NR))/2. 

ORD(NR+3)= 

. P/2 .*( l-R(NR+3)**2/P**2) 

WRITE (6,*) (R(I) , 1=1 ,NR+3) 

IF (RMAX.LE.R(NR)) WRITE (6,*) 'ERROR IN SELECTION OF RMAX ' 
IF (RMAX.LE.R(NR)) GO TO 300 
DO 6 K=1 ,NF 
6 FI (K)= 


. DCOS (DFLOAT ( 2*K - 1 ) /DFLOAT ( 2*NF )*PI / 2 . )**2 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C TIME LOOP 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

WRITE (16,*) ' irfrlrirMr Mrlririe^^ ' 


WRITE ( 16,*) ’NR,NF,NZ 3 NT0TAL,DT,P,LAM,NU,ANIS' 
WRITE (16,*) NR , NF , NZ , NTOTAL , DT , P , LAM , NU , ANI S 
WRITE (16.*) '*******■. 

WRITE(76,*) ' NR , NF , NZ , MM , DT , P , LAM , NU , ANI S ' 
WRITE (76,*) NR , NF , NZ , MM , DT , P , LAM , NU , ANI S 


— 1 — •_ .1 — 1 — > — U..I — l — U ,1 


120 


26 


364 


WR ITE (66 ,*) ' NR , NF , NZ , DT , P , LAM , NU , AN I S ' 

WRITE (66,*) NR , NF , NZ , DT , P , LAM ,NU , ANI S 
WRITE (56,*) ' NR , NF , NZ , DT , P , LAM . NU , AN I S ' 

WRITE (56,*) NR , NF , NZ , DT , P , LAM , NU , ANI S 
WR I TE ( 3 6 , * ) ' NR , NF , NZ , DT , P , LAM , NU , AN I S ' 

WRITE (36,*) NR , NF , NZ , DT , P , LAM , NU , ANI S 

WR I TE (76 "** ) * ^•’•.^>YV;;V‘VVr;V*V*V*V<V*V*VVr*ViV^Vr*VV*<VVr‘V>V‘V‘<V*V*V‘V<V*V;V‘V*V*V l :V*V/V/V‘'r‘Vi l r>wV^‘V * 

WRITE (76,*) 'PERT ON ( 1) .TIME , INT/FACE VEL. AT Rl, AT 0, SLOPE AT 
. O' 

WRITE (76 
DO 100 N=l,NTOTAL 
INDEX=0 

IF (N.GT.l) WRITE (6,*) 'WMAX' ,WMAX 
IF (N.GT.l) WRITE (16,*) 'WMAX'.WMAX 
IF (N.GT.l) WRITE (36,*) 'WMAX',WMAX 
DO 120 1=1, NR 

IF (N.EQ.l) X(I )=Z0(I )+ZDINIT(I )*DT 
IF (N.EQ.l) GO TO 120 
X(I)=Z(N-1 ,I)+ZD(N-1 ,I)*DT 
CONTINUE 

CALL ZSPOW (FCN ,NSIG , NR , ITMAX , PAR , X , FNORM ,WK, IER) 

WRITE (16 ,*) ' ************************** ' 

WRITE (16,*) 'TIME' ,N*DT 

WRITE ( 16 **’ ) * *,V V." Vr “V * V Vr ‘V * V *V Vr :V Vr Vr Vr Vr «V Vr «Y ‘.V V t V r ‘.V V? “V * V * 

WRITE ( 6 , * ) ’ ***************************** ' 

WRITE (6 ,*) 'TIME ' ,N*DT 

WRITE ( 6 , * ) ' ***************************** ’ 

WRITE (6,*) 'FNORM' , FNORM, 'INDEX' , INDEX 
WRITE (6,*)' Z0(I) 

WRITE (16,*) 'FNORM ' .FNORM , ' INDEX ’ , INDEX 
WRITE (16,*)’ Z0(I) 

WRITE (66,*) ’ ************************** 

WRITE (66 ,*) ’TIME’ , N*DT , 1 PERT ’ ,NPERT(N) 

WRITE (66 *'* ) * vfc*v*YV?#wVVwV*v*V“V#wwY“wwV#wV#V“W* * 

DO 26 1=1, NR 
Z(N, I )=X( I ) 

IF (N.GT.l) 

. ZD(N, I )=(X(I ) -Z(N-1 , I) )/DT 
IF (N.EQ.l) ZD(N,I)=ZDINIT(I) 

WRITE (66,*) R (I ) , ZD(N , I ) 

WRITE (16,*) Z0(I),X(I),X(I)-Z0(I) 

WRITE (6,*) Z0(I),X(I),X(I)-Z0(I) 

CONTINUE 

IF (NPERT(N).EQ.l) ZD(N, 1)=ZD(N, 1)+PRT 
DO 364 1=1, NR 
WRITE (6,*) R(I ) , ZD(N , I ) 

ORD(I )=X(I ) 

WR I TE (6 '** ) * ^‘wwVVwV^V^VVwV * 

ZTIPO=ZTIPN 

CALL ICSICU (R , ORD , NR+3 , BPAR , C , 30 , IER) 

VECRHO(l)=R(l) 

VECRHO (2 )=R ( 1 j+SMALL 
VECRHO ( 3 )=R ( 1 )+2*SMALL 

CALL ICSEVU (R , ORD , NR+3 , C , 30 , VECRHO , ZETA , 3 , IER) 


X(I) ’ 
X(I) ’ 



DZ 1DR i = ( ZETA ( 2 ) - ZETA ( 1 ) ) / SMALL 
DZ2DR2= (ZETA ( 3 )+ZETA ( 1 ) -2*ZETA ( 2 ) ) /SMALL**2 
ZTIPN=0RD(1 ) -R ( 1 ) "DZ 1DR1+R ( 1 )**2/2 . *DZ2DR2 
ZDO=(ZTIPN-ZTIPO)/DT 

IF (IER.GT.128) WRITE (76,*) 1 ERROR ' , IER, ' AT TIME ' ,N*DT • 

WRITE (76,760) NPERT (N ) , N*DT , ZD (N , 1 ) , ZDO , DZ 1DR 1 -R ( 1 )*DZ2DR2 
760 FORMAT (13 , 1X,F6 . 4 , IX , 3 (D16 .8 , 2X) ) 

WRITE (36,*) 'TIME ’ ,N*DT 
WRITE (56,*) 'TIME' ,N*DT 
ZERO=0 . 

DO 465 1=1 ,NR+3 
IF (I.GT.3) 

.WRITE (56,*) -R(NR+4-I), ORD ( NR+4 - I ) - Z 0 ( NR+4 - I ) 

465 WRITE (36,*) -R(NR+4-I) ,0RD(NR+4-I) 

WRITE (36,*) ZERO ,0RD(1 ) -R(1)*DZ1DR1+R( l)**2/2 .*DZ2DR2 
WRITE (56,*) ZERO , ORD ( 1 ) -R ( 1 )*DZ 1DR 1+R ( 1 )**2/ 2 . *DZ2DR2 -P/ 2 . 

WRITE (6,*) ZERO , ORD ( 1 ) -R ( 1 )*DZ 1DR 1+R ( 1 )**2/2 . *DZ2DR2 
. , ORD ( 1 ) - R ( 1 ) *DZ 1 DR 1 +R ( 1 ) ** 2 / 2 . *DZ 2DR2 - P / 2 . 

DO 466 1=1 ,NR+3 

466 WRITE (36,*) R(I),ORD(I) 

DO 467 1=1, NR 

WRITE (56,*) R ( I ) ,ORD (I ) -ZO ( I ) 

467 WRITE (6,*) R(I ) ,ORD(I ) ,ORD(I ) -20(1) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C INTERPOLATION IN CARTESIAN COORDINATES 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

332 IF (N.EQ.NTOTAL) GO TO 300 

ccccccccccccccccccccccccccccccccrcccccccccccccccccccccccccccccccccccccc 

C . CALCULATION OF U.GRAD(T) AT TIME N.DT AND AT 
C RHO=R ( J) , Z=Z (N , J) -DLOG (CSI (M) )/2 . / ( 1 . +LAM) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
IF (IVAN) GO TO 100 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C CALCULATE -(OSEEN VELOCITY) ON INTERFACE 

C (THE EXPRESSIONS FOR UR AND UZ HAVE A MINUS SIGN 

C TO SHOW THAT THEY ARE THE OPPOSITE OF THE OSEEN VEL.) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DO 820 J=1 ,NR+2 
RJ=R( J) 

ZJ=ORD( J) 

W= ( ZJ+DSQRT ( Z J**2+RJ**2 ) ) /P 
S= ( - ZJ+DSQRT ( Z J**2+R J**2 ) ) / P 
IF (S.LE.O) S=0 . 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C FLOW DOWN 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
IF (UPFLOW) GO TO 811 
UZ=- ( - 1 . +EXPIN (R2*W ) /EXPIN (R2 )+ 

. (DEXP ( -R2 ) -DEXP ( -R2*W) ) /EXPIN (R2 ) /R2/ (W+S ) ) 

UR=- ( ( 1 . -EXPIN (R2*W) /EXPIN (R2 ) )*DSQRT(S/ (S+W) ) ) 

GO TO 812 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C FLOW UP 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
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811 UZ=-(1 • -EXPIN(-R2*W)/EXPIN(R2)- 

. (DEXP ( -R2 ) -DEXP ( -R2*W ) ) /EXPIN (R2 ) /R2/ (W+S ) ) 

UR=- ( - ( 1 . -EXPIN (R2*W) /EXPIN (R2 ) )*DSQRT ( S/ (S+W) ) ) 

812 STRES ( J)=-UR 

. -DCADRE (RIGHT1 , ZEROPL , R (NR+3 ) -SMALL , 0 . , 1 . D- 7 . ERROR , IER) 

STRES (J+NR+2)=-UX 

. -DCADRE (RIGHT2 , ZEROPL , R(NR+3 ) -SMALL, 0 . , 1 . D-7 , ERROR , IER) 

DO 820 1=1 ,NR+2 

IF (I.GT.l) RD= (R (I )+R ( I - 1 ) )* . 5 
IF (I.EQ.l) RD=ZEROPL 
IF (I.LT.NR+2) RU=(R(I )+R(I+l) )*. 5 
IF (I.EQ.NR+2) RU=R (NR+3) -SMALL 
FLOW ( J , I )= 

. DCADRE (FA , RD , RU , 0 . , 1 . D ■* 7 , ERROR , IER ) 

FLOW ( J , I+NR+2 )= 

. DCADRE (FB , RD , RU , 0 . , 1 . D - 7 , ERROR , IER ) 

FLOW( J+NR+2 , I )= 

. DCADRE ( FC , RD , RU , 0 . , 1 . D - 7 , ERROR , I ER ) 

FLOW (J+NR+2, I+NR+2 )= 

. DCADRE ( FD , RD , RU , 0 . , 1 . D - 7 , ERROR , IER ) 

820 CONTINUE 

CALL LEQT1F (FLOW, 1 , 2*NR+4 , 60 , STRES ,5 ,WKAREA, IER) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CALCULATION OF INTERIOR CORRECTION TO OSEEN SOLUTION 
(WE ADD THIS SOLUTION TO THE OSEEN SOLUTION) 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
WMAX=0 . 

DO 7 J=1 , NR 
DO 7 L=1 ,NZ 
RJ=R(J) 

ZJL=Z(N,J)-DLOG(CSI(L))/2./(l.+LAM) 

ZJ=ZJL 

W= ( ZJL+DSQRT ( Z JL**2+RJ**2 ) ) /P 
WMAX=DMAX1 (W ,WMAX) 

S= ( - ZJL+DSQRT ( Z JL**2+R J**2 ) ) / P 
IF (S.LE.O) S=0. 

UR= . 5* 

. DCADRE (RIGHT1 . ZEROPL ,R (NR+3) -SMALL, 0 . , 1 . D-7 , ERROR , IER) 

UZ=.5* 

. DCADRE (RIGHT2 , ZEROPL , R (NR+3 ) -SMALL, 0 . , 1 . D- 7 , ERROR , IER) 

DO 830 1=1 ,NR+2 

IF (I.GT.l) RD= (R ( I )+R ( 1 - 1 ) )* . 5 
IF (I.EQ.l) RD=ZEROPL 
IF (I.LT.NR+2) RU=(R(I)+R(I+1) )*.5 
IF (I.EQ.NR+2) RU=R (NR+3) -SMALL 
UR=UR+.5*STRES (I )* 

. DCADRE ( FA , RD , RU , 0 . , 1 . D - 7 , ERROR , I ER ) 

.+.5*STRES( I+NR+2)* 

. DCADRE (FB. RD, RU, 0. ,1. D-7, ERROR, IER) 

UZ=UZ+ . 5*STRES ( I )* 

. DCADRE (FC . RD . RU , 0 . , 1 . D- 7 , ERROR , IER ) 

. + . 5*STRES (I+NR+2)* 

. DCADRE (FD , RD , RU , 0 . , 1 . D- 7 , ERROR , IER) 

830 CONTINUE 




ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C FLOW DOWN 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

IF (UPFLOW) GO TO 511 
WRITE (18.*) ' N= ' , N , ' Z= ' , Z JL , ' R= ' , RJ 
WRITE (18,*) ’CORRECTION' , UR, UZ 
UZ=UZ - 1 . +EXPIN (R2*V) /EXPIN (R2 )+ 

. (DEXP ( -R2 ) -DEXP ( -R2*W ) ) /EXPIN (R2 ) /R2/ (W+S ) 

UR=UR+( 1 . -EXPIN (R2*W) /EXPIN (R2 ) )*DSQRT(S/ (S+W) ) 

WRITE (18,*) ’TOTAL' , UR, UZ 
GO TO 512 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C FLOW UP 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

511 WRITE (18,*) ' N= ' , N , ' Z= ' , Z JL , ' R= ' , R J 
WRITE (18,*) ’CORRECTION’ , UR, UZ 
UZ=UZ+1 . -EXPIN (R2*W)/EXPIN(R2)- 

. (DEXP ( -R2 ) -DEXP (-R2*W)) /EXPIN (R2)/R2/ (W+S) 

UR=UR - ( 1 . -EXPIN(R2*W) /EXPIN (R2 ) )*DSQRT(S/ (S+W) ) 

WRITE (16,*) 'TOTAL' , UR, UZ 

512 DO 20 1=1, NR 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C CONTRIBUTION FROM HI & H3 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
TO=P*P I / DFLOAT ( 2*NR*NF ) *DTAN ( DFLOAT ( 2* I - 1 ) / DFLOAT ( 2*NR ) *PI / 2 . ) 
T1=T0 

.*DEXP(-Z(N , J )+Z0 (1 )+DLOG (CSI (L))/2./( 1+LAM) ) 

R00T1=DSQRT( (R(I )+R ( J) )**2+(Z (N , J) -ZO (I ) -DLOG(CSI (L) ) 

. / 2 . / ( 1+LAM ) )**2 ) 

R00T2=DSQRT ((R(I)-R(J) )**2+(Z (N , J) -ZO (I ) -DLOG ( CS I ( L ) ) 

./2./( 1+LAM) )**2) 

DO 30 K= 1 , NF 

Y=ROOT2+(ROOT1-ROOT2)*FI(K) 

Al = 

. - 2*DS I NH ( Y ) +DEXP ( Y ) *DERF ( Y/ 2 . / D S QRT ( DFLOAT ( N ) *DT ) 

, .+DSQRT(DFLOAT(N)*DT) )+ 

' . DEX P ( - Y ) *DERF ( Y / 2 . / DS QRT ( DFLOAT ( N ) *DT ) - 
. D S QRT ( DFLOAT ( N ) * DT ) ) 

A2=A1/Y 

B2= 

. - 2*DC0SH ( Y ) +DEXP ( Y ) *DERF ( Y / 2 . / DSQRT ( DFLOAT ( N ) *DT ) 

. +DS QRT ( DFLOAT ( N ) *DT ))- 
.DEXP(-Y)*DERF (Y/2 . /DSQRT (DFLOAT (N)*DT) - 
. DSQRT (DFLOAT (N)*DT)) 

.+2 . /DSQRT ( DFLOAT (N ) *DT*P I )* 

. DEXP(-Y**2/4 . /DFLOAT (N)/DT-N*DT) 

UGT(J ,L,N)=UGT(J ,L,N)+T1/DSQRT( (R00T1+Y)*(R00T2+Y) )* 

. (UZ* ( -A1+ ( B2 -A2 ) * 

. (Z(N.J)-ZO(I) -DLOG (CSI (L) ) /2 . / ( 1+LAM) )/Y) 

. +UR* (B2-A2 )*(R ( J ) - ( 0 . 5* (R00T1**2+R00T2**2 ) -Y**2 ) /2 . /R ( J) ) /Y) 

DO 40 M=1 ,NZ 

R00T3=DSQRT ( (R ( I )+R ( J) )**2+ (Z (N , J ) -ZO ( I ) 

. -DLOG (CSI (L) /CSI (H))/2./( 1+LAM) )**2) 

R00T4=DSQRT ( (R ( I ) -R ( J) )**2+ ( Z (N , J ) -ZO ( I ) 
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. -DLOG(CSI (L)/CSI(M))/2./( 1+LAM) )**2) 

0MEGA=R00T4+ (ROOTS -R00T4 )*F I (K) 

El= 

. - 2*DS INH ( OMEGA ) ■ +DEXP ( OMEGA ) *DERF ( OMEGA/ 2 . /DSQRT ( DFLOAT(N) *DT ) 

. +DSQRT (DFLOAT (N )*DT ) )+ 

. DEXP ( - OMEGA ) *DERF ( OMEGA/ 2 ./ DSQRT ( DFLOAT (N)*DT) - 
. DSQRT ( DFLOAT (N)*DT)) 

E2=E1 /OMEGA 
F2= 

. - 2*DCOSH ( OMEGA )+DEXP( OMEGA )*DERF( OMEGA/ 2 . /DSQRT (DFLOAT (N)*DT) 
.+DSQRT(DFLOAT(N)*DT) ) - 

. DEXP ( -OMEGA) *DERF (OMEGA/2 . /DSQRT (DFLOAT (N)*DT) - 
. DSQRT (DFLOAT (N )*DT) ) 

.+2 . /DSQRT ( DFLOAT (N)*DT*P I )* 

. DEXP ( -0MEGA**2/4 . /DFLOAT (N) /DT-N*DT) 

UGT(J ,L,N)=UGT(J,L,N) -TO 

. *DEXP( -Z (N, J)+ZO (I )+DLOG(CSI (L)/CSI (M) )/2 . / ( 1+LAM) ) 

. *PI*LAM/ ( 1+LAM) /DFLOAT (2*NZ ) 

.*CS2(M) 

. /DSQRT( ( ROOT3+OMEG A ) * ( R00T4+0MEG A ) ) * 

. (UZ*(-E1+(F2-E2)* 

. (Z(K.J)-ZO(I) -DLOG (CSI(L)/CSI(M))/2./( 1+LAM) ) /OMEGA) 

. +UR* (F2 -E2 )* (R ( J) - ( 0 . 5* (ROOT3**2+ROOT4**2 ) - 
. OMEGA**2 )/2 . /R ( J) ) /OMEGA) 

.*UGTO(I,H) ■ 

40 CONTINUE 
30 CONTINUE 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C CONTRIBUTION FROM H2 & H4 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO 50 NN=1,N 

IF (NN.GT.l) ZN=Z(NN-1 ,1) 

IF (NN.EQ.l) ZN=Z0 (I ) 

ROOT5=DSQRT( (R(l )+R( J) )**2+(Z(N , J)-ZN-DLOG(CSI (L) ) 

. / 2 . / ( 1+LAM) )**2 ) 

R00T6=DSQRT( (R(I )-R(J) )**2+(Z(N , J)-ZN-DLOG(CSI (L) ) 

./ 2. /(1+LAM ))** 2) 

IF (NN.EQ.l) DZDT=ZDINIT(I ) 

IF (NN.GT.l) DZDT=ZD (NN- 1 , I ) 

T2=T0/2 .*DEXP(-Z(N, J)+ZN+DLOG(CSI (L) )/2 . / (1+LAM) ) 

DO 60 K=1,NF 

Y=R00T6+ (ROOTS -ROOT6 )*FI (K) 

IF (NN.EQ.N) GO TO 68 
Cl= 

. DEXP (Y)* (DERF (Y/ 2 . /DSQRT (DFLOAT (N-NN)*DT) 

. +DSQRT ( DFLOAT ( N -NN )*DT ) ) - 

.DERF (Y/2 . /DSQRT (DFLOAT (N-NN+1)*DT) 

. +DSQRT ( DFLOAT ( N -NN+ 1 )'*DT ) ) ) + 

. DEXP ( - Y )* ( DERF ( Y/ 2 . /DSQRT ( DFLOAT ( N - NN ) *DT ) 

. -DSQRT ( DFLOAT (N-NN)*DT) ) - 
. DERF (Y/2 . /DSQRT (DFLOAT (N-NN+1 )*DT) 

. - D S QRT ( DFLO AT ( N - NN+ 1 ) * DT ) ) ) 

C2=C1/Y 


. DEXP ( Y )* ( DERF (Y/2 . / DSQRT ( DFLOAT ( N - NN ) *DT ) 

. +DSQRT ( DFLOAT ( N - NN ) *DT) ) - 
. DERF ( Y/ 2 . /DSQRT ( DFLOAT (N-NN+1) *DT ) 

. +DSQRT (DFLOAT (N-NN+1 )*DT) ) ) - 
. DEXP ( - V DERF ( V /:./ DSQRT ( DFLOAT ( N -NN ) *DT ) 

. -DSQRT (DFLOAT(N-NN)*DT) ) - 
. DERF ( Y / 2 . / DSQRT ( DFLOAT ( N-NN+ 1 ) *DT ) 

. - DSQRT ( DFLOAT ( N - NN+ 1 ) * DT ) ) ) + 

. 2 . /DSQRT (PI*DFLOAT(N-NN)*DT) 

. *DEXP ( - Y** 2 / 4 . / DFLO AT ( N - NN ) / DT - ( N - NN ) *DT ) - 
. 2 . /DSQRT ( P I DFLOAT ( N -NN+ 1 ) *DT ) 

. *DEXP ( - Y**2 / 4 . / DFLOAT ( N -NN+ 1) / DT - ( N -NN+ 1 ) *DT ) 

UGT (J,L, N)=UGT(J , L , N ) +T2* ( 2 . +DZDT ) / 

. DSQRT ( ( ROOTS +Y ) * ( ROOT 6+Y ) ) * 

. (UZ*(-C1+(D2-C2)* 

. (Z(N, J)-ZN-DLOG(CSI(L))/2./( 1+LAM) )/Y) 

. +UR*(D2-C2)*(R ( J ) - (0 . 5* (ROOT5**2+ROOT6**2 ) -Y**2 ) /2 . /R( J) )/Y) 
GO TO 69 

68 Cl= 

.DEXP(Y)* 

.DERFC(Y/2 . /DSQRT(DT) 

.+DSQRT(DT) )+ 

. DEXP (-Y )••'•• 

. DERFC (Y/2 . /DSQRT (DT) 

. -DSQRT(DT) ) 

C2=C1/Y 

D2= 

. DEXP (Y )•■'•• 

. DERFC (Y/2. /DSQRT(DT) 

,+DSQRT(DT) ) - 
. DEXP (-Y )-•■•- 

. DERFC (Y/2. /DSQRT(DT) 

. -DSQRT(DT))- 
. 2 . /DSQRT (PI*DT) 

. •••DEXP ( - Y----2/4 . /DT-DT) 

UGT ( J , L , N ) =UGT ( J , L , N ) +T2* ( 2 . +DZDT ) / 

. DSQRT ( (ROOT5+Y)*(ROOT6+Y) )* 

. (UZ*(-C1+(D2-C2)* 

. (Z (N , J ) -ZN-DLOG (CS1 (L) )/2 . / ( 1+LAM) ) /Y) 

.+UR*(D2-C2)*(R(J)-(0 . 5*(ROOT5**2+ROOT6**2) -Y**2)/2 . /R(J) )/Y) 

69 DO 70 M=1 ,NZ 

IF (NN.EQ.l) UGTOLD=UGTO ( I . M ) 

IF (NN.GT.l) UGTOLD=UGT ( I , M , NN - 1 ) 

R00T7=DSQRT ( (R ( I )+R ( J ) )**2+ ( Z (N , J ) - ZN 
.+DLOG(CSI (M)/CSI (L) )/2 . / (1+LAM) )**2) 

R00T8=DSQRT ( (R( I ) -R( J) )**2+(Z (N , J) -ZN 
.+DLOG(CSI (M)/CSI (L) ) / 2 . / (1+LAM) )**2) 

0MEGA=R00T8+ (R00T7 -R00T8 )*FI (K) 

IF (NN.EQ.N) GO TO 71 
Ql= 

. DEXP ( OMEGA ) ■ * ( DERF ( OMEGA / 2 . / DSQRT ( DFLOAT ( N - NN ) *DT ) 

. +DSQRT (DFLOAT (N -NN ) '•'••DT) ) - 
.DERF(0MEGA/2 . /DSQRT (DFLOAT (N-NN+1 )*DT) 

. +DSQRT(DFLOAT (N-NN+1 )*DT) ) )+ 


. DEXP ( -OMEGA)* (DERF (OMEGA/2 . /DSQRT(DFLOAT(N-NN)*DT) 

. -DSQRT(DFLOAT(N-NN)*DT) ) - 
. DERF ( OMEGA/ 2 ./ DSQRT ( DFLOAT (N-NN+1 ) *DT ) 

. -DSQRT (DFLOAT (N-NN+ 1 )*DT) ) ) 

Q2=Ql/OMEGA 

S2= 

. DEXP (OMEGA)* (DERF (OMEGA/ 2 . /DSQRT (DFLOAT (N-NN)*DT) 

. +DSQRT ( DFLOAT (N-NN)*DT) ) - 

.DERF (OMEGA/ 2 . / DSQRT ( DFLOAT ( N -NN+ 1 )*DT) 

. +DSQRT (DFLOAT (N-NN+1 )*DT) ))- 

. DEXP (- OMEGA )*( DERF ( OMEGA / 2 ./ DSQRT ( DFLOAT ( N -NN ) *DT ) 

. -DSQRT (DFLOAT (N-NN)*DT) ) - 

.DERF (OMEGA/2 . /DSQRT ( DFLOAT ( N - NN+ 1 ) *DT ) 

. -DSQRT(DFLOAT(N-NN+l )*DT) ) )+ 

. 2 . / DSQRT ( P I*DFLOAT ( N -NN ) *DT ) 

. *DEXP ( - OMEG A** 2 / 4 . / DF LO AT ( N - NN ) / DT - ( N - NN ) *DT ) - 
. 2 . / D S QRT ( P I *DF LOAT (N-NN+1 ) *DT ) 

. *DEXP ( -0MEGA**2/4 . /DFLOAT (N-NN+1 ) /DT- (N-NN+1 )*DT) 
UGT(J,L,N)=UGT(J,L,N) 

. -T2*PI*LAM/ ( 1+LAM) /DFLOAT (NZ)*CS I (M)** (-.5/(1. +LAM ) ) 
•*CS2 (M) 

. /DSQRT ( (ROOT7+OMEGA)* (R00T8+0MEGA) )* 

. (UZ*(-Q1+(S2-Q2)* 

. (Z(N,J)-ZN-DLOG(CSI(L)/CSI(M))/2./( 1+LAM)) /OMEGA) 

. +UR"'(S2 -Q2 )••• 

. (R ( J) - ( 0 . 5* (ROOT7** 2+ROOT8** 2 ) -OMEGA**2 ) /2 . /R ( J) ) /OMEGA) 
. *UGTOLD 
GO TO 70 
Ql= 

.DEXP (OMEGA)* 

.DERFC (OMEGA/ 2. /DSQRT (DT) 

.+DSQRT(DT) )+ 

.DEXP (-OMEGA)* 

. DERFC (OMEGA/2 . /DSQRT (DT) 

. -DSQRT (DT)) 

Q2=Q1 /OMEGA 
S2= 

. DEXP (OMEG A)* 

. DERFC (OMEGA/ 2 . /DSQRT (DT) 

.+DSQRT(DT) ) - 
•DEXP (-OMEGA)* 

.DERFC (OMEGA/ 2 . /DSQRT (DT) 

. -DSQRT(DT))- 
. 2 . /DSQRT(PI*DT) 

. *DEXP ( -0MEGA**2/4 . /DT-DT) 

UGT(J ,L,N)=UGT(J,L,N) 

. -T2*PI*LAM/ ( 1+LAM ) /DFLOAT (NZ )*CSI (M)**(- . 5/ (1 .+LAM) ) 
•*CS2(M) 

. /DSQRT ( (ROOT7+OMEGA)*(ROOT8+OMEGA) )* 

. (UZ*(-Q1+(S2-Q2)* 

. (Z (N, J) -ZN-DLOG(CSI (L)/CSI(M))/2. /( 1+LAM) ) /OMEGA) 

. +UR*(S2-Q2)* 

. (R(J)-(O. 5* (ROOT7**2+ROOT8**2 ) -OMEGA**2 ) / 2 . /R(J) ) /OMEGA) 
,*UGTOLD 


70 CONTINUE 

60 CONTINUE OR 1GM&H E©3g jg 

50 CONTINUE OF POOR QUALITY! 

20 CONTINUE 1 

7 CONTINUE 
100 CONTINUE 
300 STOP 
END 

SUBROUTINE FCN(X,F,NR.PAR) 

IMPLICIT REAL- 8 (A-H.0-Z) 

REAL*8 X(NR),F(NR),PAR(1), 

. UGT(30 , 30, 1000 ) ,UGT0 (30 , 30) , Z ( 1000 ,30), 

. ZD( 1000 ,30) ,Z0(30) ,R(30) ,CSI (30) ,CS2 (30) ,FI (30) , ZDINIT(30) , 
.0RD(30) ,VECRHO(3) ,ZETA(3) ,BPAR(4) ,C(30,3) ,LAM,NU 
LOGICAL IVAN 

COMMON /A/UGT , UGTO , Z , ZD , ZDINIT , CSI , CS2 , FI , ZO , LAM , PI , 

. DT , NU , AN I S , NF , N Z , NT , I NDEX , I VAN 
. /B/C ,R ,0RD ,BPAR 
./D/ SMALL, P 
F ( NR )=X (NR ) - ZO (NR ) 

I NDEX= I NDEX+ 1 
IF (NU.EQ.O.) GO TO 999 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C CONTRIBUTION FROM CURVATURE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
0RD(NR+1 )= 

. P/2 . "•( 1 -R(NR+1 ) ** 2 / P" " 2 ) 

ORD(NR+2)= 

. P / 2 . * ( 1 - R (NR+2 ) * 2 / P- • " 2 ) 

0RD(NR+3)= 

. P/2 . * ( 1 -R(NR+3)**2/P**2) 

DO 4 1=1. NR 
4 ORD(I )=X( I ) 

CALL . I C S I CU ( R , ORD , NR+3 . BPAR , C , 30 , IER ) 

DO 992 J=1 ,NR- 1 
RJ=R( J) 

IF (J.GT. 1) GO TO 998 
VECRHO ( 1 )=RJ 
VECRHO ( 2 )=RJ+SMALL 
VECRHO ( 3 ) =R J+SMALL* 2 

CALL ICSEVU (R . ORD , NR+3 . C . 30 , VECRHO , ZETA , 3 , IER ) 

DZ1DR1=(ZETA (2)-ZETA( 1) ) /SMALL 
DZ2DR2=(ZETA(3)+ZETA( 1 ) -2*ZETA (2) )/SMALL**2 
GO TO 997 

998 VECRHO (1)=RJ -SMALL 
VECRHO ( 2 )=RJ 
VECRHO ( 3 )=R J+SMALL 

CALL I CSEVU (R , ORD , NR+3 , C , 30 , VECRHO , ZETA , 3 , IER ) 

DZ1DR1=(ZETA (3) -ZETA ( 1 ) )/2 . /SMALL 
DZ2DR2=( ZETA ( 3 ) +ZETA ( 1 ) - 2*ZETA ( 2 ) ) / SMALL**2 
997 F(J)=-NU*(RJ*DZ2DR2+DZ1DR1**3+DZ1DR1)/RJ 
. / ( 1+DZ 1DR 1 "••'•'2 ) - '•• 1 . 5 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C ANISOTROPY 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

.*( 1 -ANI S*8*DZ 1DR 1**2 / 

. (1+DZ1DR1**2)**2) 

992 CONTINUE 
999 DO 2 J=1 ,NR- 1 

IF (NU.EQ.O.) F ( J )=0 . 

DO 20 1=1, NR 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C CONTRIBUTION FROM -DEL 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

RPLUS=DSQRT( 

. (R(I )+RJ)**2+(Z0( J)-ZO(I ) )**2) 

RMIN=DSQRT( 

. (R ( I ) -RJf )**2+ ( ZO ( J) -ZO ( I ) )**2 ) 

DO 5 K=1 ,NF 

Y=RMIN+(RPLUS-RMIN)*FI (K) 

F ( J ) =F ( J ) - P* P I / DFLO AT ( NR*NF ) *DTAN 
. ( DFLOAT ( 2* I r 1 ) /DFLOAT ( 2*NR ) *PI / 2 . ) 

-*DEXP.(ZO(I)-ZO(J)-Y)/ 

. DSQRT ( ( RM I N+Y ) * ( RPLUS+Y ) ) 

5 CONTINUE 

DO 10 M=1 ,NZ 

R E S 1 = P I ** 2* P / 2 . /DFLOAT (NR*NF*NZ)*LAM/ ( 1 ,+LAM)* 

. DTAN ( DFLOAT (2*1-1)/ DFLOAT ( 2*NR ) *P I / 2 . ) 

.*CS2 (M) 

.*CSI (M )**( - 1 . /2 . / ( 1 . +LAM ) ) 

.*DEXP(ZO (I ) -ZO ( J) ) 

,*UGTO(I,M) 

RPLUS=DSQRT ( (R ( I )+RJ )**2+ (ZO ( J ) -ZO ( I ) 

. +DLOG (CSI (M) )/2 . / ( 1+LAM) ).**2 ) 

RM I N=DSQRT ( ( R ( I ) - R J ) ** 2+ ( Z 0 ( J ) - Z 0 ( I ) 

.+DLOG (CSI (M) )/2 . / ( 1 +LAM ) ) ** 2 ) 

DO 10 K= 1 , NF 

OMEG A=RM IN+ (RPLUS - RM IN ) *F I ( K ) 

F(J)=F(J)+RES1*DEXP( -OMEGA)/ 

. DSQRT ( (RPLUS+OMEGA)*(RMIN+OMEGA) ) 

10 CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C CONTRIBUTION FROM II & 13 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
TO=P*P 1 / DFLOAT ( 2*NR*NF ) * DTAN ( DFLOAT (2*1-1) /DFLOAT ( 2*NR ) *P 1/2 . ) 
T1=T0 

.*DEXP( -X( J)+ZO ( I ) ) 

ROOT 1=DSQRT ( ( R ( I ) +R J )**2+ ( X ( J ) - ZO ( I ) ) **2 ) 

R00T2=DSQRT ( (R ( I ) -RJ )**2+ (X ( J ) -ZO ( I ) )**2 ) 

DO 30 K= 1 , NF 

Y=R00T2+ ( ROOT 1 - R00T2 ) *F I ( K ) 

F ( J)=F (J)+T1 /DSQRT ( (R00T1+Y)*(R00T2+Y) )* 

. ( - 2*DS INH ( Y ) +DEXP ( Y ) *DERF ( Y / 2 . /DSQRT ( DFLOAT ( NT ) *DT ) 

. +DSQRT ( DFLOAT ( NT ) *DT ) )+ 

.DEXP(-Y)*DERF(Y/2 . / DSQRT ( DFLOAT ( NT )*DT) - 
. DSQRT (DFLOAT ( NT ) *DT ) ) ) 

IF (IVAN) GO TO 30 
DO 40 M=1,NZ 


•V 


R00T3=DSQRT ( (R ( I )+RJ)**2+(X ( J) -ZO ( I ) 

. +DLOG ( CS 1 (M ) ) / 2 . / ( 1 +LAM ) ) ** 2 ) 

R00T4=DSQRT ( ( R ( I ) - R J ) ** 2+ ( X ( J ) - Z 0 ( I ) 

. +DLOG (CS 1 (M ) ) / 2 . / ( 1+LAM ) ) **2 ) 

0MEGA=R00T4+ (ROOTS - R00T4 )*FI (K) 

F ( J ) =F ( J ) -T 1* P I "LAM / (1+LAM) /DFLOAT ( 2*KZ ) •'• CS I (M) ** (- . 5/ (1 .+LAM) ) 
.*CS2 (M) 

. /DSQRT ( (R00T3+0MEGA )■•'•' (ROOT4+OMEGA ) )* 

. ( - 2*DS INH ( OMEGA ) +DEXP ( OMEGA )*DERF (OMEGA/ 2 . / DSQRT ( DFLOAT (NT ) *DT ) 

. +DSQRT ( DFLOAT ( NT )*DT) )+ 

. DEXP ( - OMEGA ) *DERF ( OMEGA/ 2 . / DSQRT (DFLOAT ( NT ) *DT ) - 
. DSQRT (DFLOAT (NT )*DT) ) ) 

.*UGTO(I,M) 

40 CONTINUE 
30 CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C CONTRIBUTION FROM 12 & 14 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DO 50 N= 1 , NT 
IF (N.EQ.NT) ZN=X(I ) 

IF (N.LT.NT) ZN=Z(N,I) 

R00T5=DSQRT ( (R ( I ) +R J )**2+ (X ( J ) - ZN ) **2 ) 

R00T6=DSQRT ( (R ( I ) -RJ)**2+(X( J) -ZN)**2) 

IF (N . EQ . NT . OR . NT . EQ . 1 ) GO TO 57 
DZDT=ZD(N, I ) 

GO TO 56 

57 IF (NT.EQ.l) DZDT=ZDINIT(I ) 

IF (NT.GT.l) DZDT=(X(I)-Z(N-1,I))/DT 
56 T2=TO/2.*DEXP(-X(J)+ZN) 

DO 60 K= 1 , NF 

Y=ROOT6+ ( ROOT5 - ROOT 6 ) *F I (K ) 

IF (N.NE.NT) 

. F ( J )=F ( J )+T2*(2 . +DZDT) / 

. DSQRT ( ( ROOT5+ Y ) * (ROOT 6+Y ) )* 

. ( DEXP ( Y )-•'•- (DERF(Y/ 2. /DSQRT (DFLOAT (NT-N)*DT) 

. +D S QRT ( DFLO AT ( NT - N ) * DT ) ) - 
. DERF ( Y / 2 . /DSQRT ( DFLOAT ( NT - N+ 1 ) *DT ) 

. +DSQRT(DFLOAT(NT-N+l )*DT) ) )+ 

. DEXP ( - Y )•■'••( DERF ( Y / 2 . / DS QRT ( DFLOAT ( NT - N ) •■'••DT ) 

. -DSQRT(DFLOAT(NT-N)*DT) ) - 
. DERF ( Y / 2 . /DSQRT ( DFLOAT ( NT - N+ 1 ) * DT ) 

. -DSQRT(DFLOAT(NT-N+l )*DT) ) ) ) 

IF (N.EQ.NT) 

. F(J)=F(J)+T2* (2 .+DZDT)/ 

.DSQRT( (ROOT5+Y)* (ROOT6+Y) )* 

. (DEXP(Y)* 

. DERFC (Y/ 2 . /DSQRT (DT) 

.+DSQRT(DT))+ 

.DEXP(-Y)* 

. DERFC (Y/2. /DSQRT (DT) 

. -DSQRT (DT) ) ) 

IF (IVAN) GO TO 60 
IF (N.EQ.l) ZOLD=ZO(I ) 

IF (N.GT.l) ZOLD=Z (N-l , I ) 
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IF (N.EQ.NT) GO TO 78 
DO 70 M=1 ,NZ 

IF (N.EQ.l) UGTOLD=UGTO ( I , M) 

IF (N.GT.l) UGTOLD=UGT (I , M , N- 1 ) 

ROOT 7=DSQRT ( ( R ( I ) +R J ) **2+ ( X ( J > r ZOLD 
. +DLOG ( CS I (M) ) / 2 . / ( 1+LAM ) )**2 ) 

R00T8=DSQRT( (R ( I ) -RJ)**2+(X(J) -ZOLD 
.+DLOG(CSI (M) )/2 . / ( 1+LAM) )**2) 

OMEG A=R00T8+ ( ROOT? -ROOTS )*FI (K ) 

F ( J)=F ( J ) - 

.T0/2.*DEXP(-X(J)+Z0LD) 

. *PI*LAM/ ( 1+LAM) /DFLOAT(NZ)*CSI (M)**(- . 5/(1. +LAM ) ) 

. *CS2(M) 

. / DSQRT ( (R00T7 +OMEG A ) * ( R00T8+0MEG A ) )* 

. (DEXP(0MEGA)*(DERF(0MEGA/2./DSQRT(DFL0AT(NT-N)*DT) 

. +DSQRT ( DFLOAT ( NT- N )*DT) ) - 
. DERF ( OMEG A/ 2 . / DSQRT ( DFLOAT ( NT - N+ 1 ) *DT ) 

. +DSQRT ( DFLOAT (NT-N+1 )*DT) ) )+ 

.DEXP( -OMEGA)* (DERF (OMEGA/2 . /DSQRT(DFLOAT(NT-N)*DT) 

. -DSQRT (DFLOAT (NT-N)*DT) ) - 
. DERF ( OMEGA / 2 . /DSQRT ( DFLOAT ( NT - N+ 1 ) *DT ) 

. -DSQRT (DFLOAT (NT-N+1 )*DT) ) ) ) 

. '-UGTOLD 
70 CONTINUE 
GO TO 60 

78 DO 75 M=1 ,NZ 

IF (N.EQ.l) UGT0LD=1 ! GT0 (I . M ) 

IF (N.GT.l) UGTOLD=UGT ( I , M , N - 1 ) 

R00T7=DSQRT ( (R ( I )+R J )**2+ (X ( J ) - ZOLD 
. +D LOG ( C S I ( M ) ) / 2 . / ( 1 + LAM ) ) ** 2 ) 

R00T8=DSQRT ( ( R ( I ) - R J ) **2+ ( X ( J ) - ZOLD 
. +DLOG ( CS I ( M ) ) / 2 . / ( 1 +LAM ) )** 2 ) 

0MEGA=R00T8+ (ROOT? -ROOTS )*F1(K) 

' F (J)=F(J) - 

. TO / 2 . *DEXP ( -X ( J ) +ZOLD ) 

. *PI*LAM/ ( 1 . +LAM ) / DFLOAT ( N Z ) *C S I (M )** (-.5/(1. +LAM) ) 

.*CS2 (M) 

. /DSQRT ( (R00T7+0MEGA )* (R00T8+0MEGA ) )* 

. (DEXP (OMEGA)* 

. DERFC ( OMEGA/ 2 . /DSQRT (DT) 

.+DSQRT(DT) )+ 

.DEXP (-OMEGA)* 

.DERFC (OMEGA/ 2. /DSQRT (DT) 

. -DSQRT(DT) ) ) 

. *UGTOLD 
75 CONTINUE 
60 CONTINUE 
50 CONTINUE 
20 CONTINUE 
2 CONTINUE 
RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
FUNCTION EXPIN(P) 




# 


139 


OSIGMAL' page b 

°E POOR QUALITY 


f 


IMPLICIT REALMS (A-H.O-Z) 

IF (P.GT. 14. ) GO TO 30 

EXPIN=-0 .57721566490245707 l-DLOG(P) 

DOE=-l . 

DO 20 1=1.1000000 
DOE=-DOE*P/DFLOAT(I ) 

IF (DABS (DOE) /DFLOAT ( I ) . LE . 1 . D- 16)G0 TO 40 
20 EXPIN=EXPIN+DOE/DFLOAT( 1 ) 

30 EXPIN=DEXP( -P-DLOG(P) )*( 1 - 1 . /P+2/P**2-6/P**3 
,+24./P**4) 

40 RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

FUNCTION EXPIN2 (P) 

IMPLICIT REALMS (A-H.O-Z) 

EXTERNAL EXPIN 
EXPIN2=DEXP( -P) -P*EXPIN (P) 

RETURN 

END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

FUNCTION FA(RHO) 

IMPLICIT REAL*S (A-H.O-Z) 

REAL*S QUAD(30) ,VECRH0(3) ,ZETA(3) ,C(30,3) ,R(30) ,ORD(30) 

COMMON 

./B/C.R.ORD.BPAR 
. /D/SMALL. P 

. /C/QUAD. RJ.ZJ.R2. MM. NR 
IF (RHO.GT. (R( 1 )+SMALL) ) GO TO 996 
VECRHO( 1 )=R( 1 ) 

VECRHO ( 2 )=R ( 1)+SMALL 
VECRH0(3)=R( 1 )+SMALL*2 

CALL ICSEVU(R, ORD. NR+3.C, 30, VECRHO, ZETA, 3, IER) 

DER 1= ( ZETA ( 2 ) - ZETA ( 1 ) ) / SMALL 
DER2= ( ZETA ( 3 ) +ZETA ( 1 ) - 2*ZETA ( 2 ) ) / SMALL** 2 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C ZETA (2) WILL NOW REPRESENT ZETA IN [ ZEROPL.R( 1) ] 

C DZ1DR1 WILL NOW REPRESENT DER IV IN [ZEROPL.R(l) ] 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
ZETA ( 2 ) =ZETA ( 1 ) + ( RHO - R ( 1 ) ) *DER 1 + ( RHO - R ( 1 ) ) **2 / 2 . *DER2 
DZ1DR1=DER1+(RH0-R( 1 ) )*DER2 • 

GO TO 997 

998 VECRHO (1)=RH0- SMALL 
VECRHO ( 2 )=RHO 
VECRHO (3 )=RHO+SMALL 

CALL ICSEVU(R, ORD. NR+3.C, 30. VECRHO. ZETA, 3, IER) 
DZ1DR1=(ZETA(3)»ZETA(1) )/2. /SMALL 
997 B=2*RH0*RJ 

A=RH0**2+R J**2+ ( ZETA ( 2 ) - Z J ) **2 
FA=0 . 

DO 10 M=1 ,MM 

DIST=QUAD(M)*(DSQRT (A+B) -DSQRT(A-B) )* . 5 
.+(DSQRT (A+B)+DSQRT (A-B) )* . 5 
THC0S=(A-DIST**2)/B 

10 FA=FA+1 . /DFLOAT (MM)*RHO*DSQRT ( 1 . +DZ 1DR 1**2 ) 
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./DSQRT((DSQRT(A-B)+DIST)*(DSQRT(B+A)+DIST)) 

. * (THCOS+ (R J-RHO*THCOS )* (RJ*THCOS-RHO)/DI ST**2 ) 

RETURN . 

END 

FUNCTION FB(RHO) 

IMPLICIT REALMS (A-H.O-Z) 

REAL* 8 QUAD (30) ,VECRHO(3) .ZETA (3) ,C(30,3) ,R(30) ,ORD(30) 

COMMON 

./b/c,r,ord.bpar 

./D/ SMALL, P 

. /C/QUAD ,RJ , ZJ ,R2 , MM , NR 
IF (RHO.GT. (R ( 1 )+SMALL) ) GO TO 998 
VECRHO( 1 )=R ( 1 ) 

VECRHO ( 2 )=R ( 1 ) +SMALL • 

VECRH0(3)=R( 1 )+SMALL*2 

CALL ICSEVU (R , ORD , NR+3 , C , 30 , VECRHO , ZETA , 3 , IER) 

DER 1 = ( ZETA ( 2 ) - ZETA ( 1 ) ) / SMALL 
DER2= ( ZETA ( 3 ) +ZETA ( 1 ) - 2*ZETA ( 2 ) ) / SMALL**2 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
ZETA(2) WILL NOW REPRESENT ZETA IN [ ZEROPL,R( 1 ) ] 

DZ1DR1 WILL NOW REPRESENT DER IV IN ( ZER0PL,R(1) ] 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
ZETA ( 2 ) =ZETA ( 1 ) + ( RHO - R ( 1 ) ) *DER 1+ ( RHO -R ( 1 ) ) **2 / 2 . *DER2 
DZ1DR1=DER1+(RH0-R ( 1 ) )*DER2 
GO TO 997 

998 VECRHO ( 1 )=RHO- SMALL 
VECRHO ( 2 )=RHO 
VECRHO ( 3 )=RHO+SMALL 

CALL ICSEVU (R , ORE , NR+3 , C , 30 , VECRHO , ZETA , 3 , IER) 

DZ1DR1=(ZETA(3) -ZETA ( 1 ) ) / 2 . /SMALL 
99 7 B=2*RHO*R.l 

A=RHO**2+R J**2+ ( ZETA ( 2 ) - Z J )**2 
FB=0 . 

DO 10 M= 1 , MM 
DIST=QUAD(M) 

. * ( DSQRT ( A+B ) -DSQRT ( A - B ) ) * . 5 
.+( DSQRT ( A+B ) +DSQRT ( A - B ) ) * . 5 
THCOS= (A -DI ST**2 ) / B 

10 FB=FB+1 . /DFLOAT (MM ) *RHO*DSQRT ( 1 . +DZ 1DR 1**2 ) 

. /DSQRT( (DSQRT(A-B)+DIST)*(DSQRT(B+A)+DIST) ) 

. * (RJ-RHO*THCOS ) /DI ST* ( ZJ -ZETA (2 ) ) /DIST 
RETURN 
END 

FUNCTION FC (RHO) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL* 8 QUAD (30) .VECRHO (3) ,ZETA(3) ,C(30,3) ,R(30) ,ORD(30) 

COMMON 

. /B/C .R.ORD.BPAR 
. /D/SMALL, P 

. /C/QUAD, RJ.ZJ.R2, MM, NR 
IF (RHO.GT. (R(1)+SMALL)) GO TO 998 
VECRHO ( 1 )=R( 1 ) 

VECRHO ( 2 )=R ( 1 ) +SMALL 
VECRHO ( 3 )=R ( 1 )+SMALL*2 
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CALL I CSEVU (R , ORD , NR+3 , C , 30 , VECRHO , ZETA , 3 , IER) 

DER 1= ( ZETA ( 2 ) -ZETA ( 1 ) ) /SMALL 

DER2= ( ZETA ( 3 ) +ZETA ( 1 ) - 2*ZETA ( 2 ) ) /SMALL**2 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C ZETA (2) WILL NOW REPRESENT ZETA IN [ ZEROPL ,R ( 1 ) ] 

C DZ1DR1 WILL NOW REPRESENT DER IV IN (ZEROPL,R(l) ] 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ZETA ( 2 )=ZETA ( 1 )+ (RHO-R ( 1 ) )*DER 1+ ( RHO-R ( 1 ) )**2 / 2 . *DER2 
DZ 1DR 1=DER 1+ ( RHO-R ( 1 ) )*DER2 
GO TO 997 

998 VECRHO (l)=RHO-SMALL 
VECRHO ( 2 )=RHO 
VECRHO ( 3 )=RHO+SMALL 

CALL I CSEVU (R , ORD , NR+3 , C , 30 , VECRHO , ZETA , 3 , IER ) 

DZ1DR1=(ZETA( 3) -ZETA( 1 ) ) / 2 . /SMALL 

997 B=2*RHO*RJ 

A=RHO**2+RJ**2+ ( ZETA (2 ) - Z J )**2 
FC=0 . 

DO 10 M=1 ,MM 
DIST=QUAD(M) 

. * ( DSQRT ( A+B ) - DSQRT ( A -B ) ) * . 5 
.+(DSQRT (A+B)+DSQRT(A-B) )* . 5 
THCOS=(A-DIST**2)/B 

10 FC=FC+ 1 . / DFLOAT ( MM ) *RHO*DSQRT (1 .+DZ1DR 1**2 ) 

. /DSQRT ( (DSQRT (A-B)+DIST)* (DSQRT(B+A)+DIST) ) 

. * ( RJ*THCOS -RHO ) / D I ST* ( Z J - ZETA ( 2 ) ) / D I ST 
RETURN 
END 

FUNCTION FD 'RHO) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 QUAD (30) , VECRHO (3) , ZETA (3) ,C(30 , 3) ,R(30) ,ORD(30) 

COMMON 

. / B / C , R , ORD , BPAR 
./D/ SMALL, P 

. /C/QUAD ,RJ ,ZJ ,R2 ,MM,NR 
IF (RHO.GT. (R(1)+SMALL)) GO TO 998 
VECRHO ( 1)=R( 1 ) 

VECRHO ( 2 )=R(1)+SMALL 
VECRHO ( 3 )=R(1)+SMALL*2 

CALL ICSEVU (R . ORD , NR+3 , C , 30 , VECRHO , ZETA , 3 , IER ) 

DER 1 = ( ZETA ( 2 ) - ZETA ( 1 ) ) / SMALL 
DER2= ( ZETA ( 3 )+ZETA ( 1 ) - 2*ZETA (2 ) )/SMALL**2 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C ZETA (2) WILL NOW REPRESENT ZETA IN [ ZEROPL,R( 1 ) ] 

C DZ1DR1 WILL NOW REPRESENT DERIV IN [ZEROPL,R( 1 ) ] 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
ZETA ( 2 ) = ZET A ( 1 ) + ( RHO - R ( 1 ) ) *DER 1+ (RHO - R ( 1 ) ) •“'• 2 / 2 . *DER 2 
DZ 1DR 1=DER 1+ (RHO-R ( 1 ) )*DER2 
GO TO 997 

998 VECRHO (l)=RHO- SMALL 
VECRHO (2 )=RHO 
VECRHO ( 3 )=RHO+SMALL 

CALL I CSEVU (R , ORD , NR+3 , C , 30 , VECRHO , ZETA , 3 , I ER ) 

DZ 1DR 1= ( ZETA ( 3 ) -ZETA ( 1 ) ) / 2 . / SMALL 
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B=2*RH0*RJ 

A=RH0**2+R J**2+ ( ZETA ( 2 ) -Z J )**2 
FD=0 . 

DO 10 M= 1 , MM 
DIST=QUAD(M) 

. * ( DSQRT ( A+B ) -DSQRT ( A - B ) )* . 5 
.+ (DSQRT (A+B)+DSQRT(A-B) )*. 5 
THCOS= ( A - D I ST* "2 ) / B 
10 FD=FD+1 . /DFLOAT ( MM ) *RHO*DSQRT ( 1 . +DZ 1DR 1**2 ) 

. / DSQRT ( ( DSQRT ( A - B ) +D I ST ) * ( DSQRT ( B+A ) +D I ST ) ) 

. * ( 1+ ( ( ZETA (2 ) -Z J ) /DI ST)**2 ) 

RETURN’ 

END 

FUNCTION RIGHT1 (RHO) 

IMPLICIT REAL*8 (A-H.O-Z) 

REAL*8 QUAD (30) ,VECRH0(3) ,ZETA(3) ,C(30 , 3) ,R(30) ,ORD(30) 

LOGICAL UPFLOW 
COMMON 

. /B/C , R , ORD , BPAR 
./D/ SMALL, P 

. /C/QUAD , RJ , ZJ , R2 , MM , NR , UPFLOW 
EXTERNAL EXPIN 

IF (RHO.GT. (R(1)+SMALL)) GO TO 998 
VECRHO( 1 )=R( 1 ) 

VECRHO (2 )=R ( 1 )+SMALL 
VECRH0(3)=R( 1)+SMALL*2 

CALL ICSEVU ( R , ORD . NR+3 , C , 30 , VECRHO , ZETA . 3 , IER) 

DER 1= ( ZETA ( 2 )- ZETA ( 1 ))/ SMALL 
DER2= ( ZETA ( 3 )■ +ZETA ( 1 j - 2*ZETA ( 2 ) ) / SMALL**2 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C ZETA (2) WILL NOW REPRESENT ZETA IN [ZEROPL,R(l) ] 

C DZ1DR1 WILL NOW REPRESENT DER IV IN ( ZEROPL,R( 1) ] 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ZETA ( 2 ) =ZETA ( 1 ) + ( RHO - R ( 1 ) ) *DER 1+ ( RHO -R ( 1 ) ) **2 / 2 . *DER2 
DZ1DR1=DER1+(RH0-R( 1) )*DER2 
GO TO 997 

998 VECRHO( 1 )=RHO-SMALL 
VECRHO ( 2 )=RHO 
VECRHO ( 3 )=RHO+SMALL 

CALL ICSEVU (R . ORD , NR+3 , C , 30 , VECRHO , ZETA , 3 , IER ) 

DZ 1DR 1 = ( ZETA ( 3 ) - ZETA ( 1 ) ) / 2 . / SMALL 
997 W= ( ZETA ( 2 ) +DSQRT ( ZETA ( 2 ) **2+RHO**2 ) ) / P 
S= ( - ZETA ( 2 ) +DSQRT ( ZETA ( 2 )**2+RHO**2 ) ) / P 
IF (S.LE.O) S=0. 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C FLOW DOWN 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
IF (UPFLOW) GO TO 811 
UZ=- (-1 . +EXPIN (R2*W) /EXPIN (R2)+ 

. (DEXP ( -R2 ) -DEXP ( -R2*W ) ) /EXPIN (R2 ) /R2/ (W+S ) ) 

UR=- ( ( 1 . -EXPIN (R2*W) /EXPIN(R2 ) )*DSQRT(S/ (S+W) ) ) 

GO TO 812 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C FLOW UP 
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Sii UZ=-(1 . -EXPIN (R2*W) /EXPIN (R2 ) - 

. (DEXP(-R2) -DEXP( -R2*W) ) /EXPIN (R2)/R2/ (W+S) ) 

UR=- ( - ( 1 . -EXPIN (R2*W ) / EXP I N ( R2 ) ) *DSQRT ( S / ( S+W ) ) ) 

812 B=2*RH0*RJ 

A=RH 0 2 +RJ 2+ ( ZETA ( 2 ) - Z J )**2 
RIGHT1=0 . 

DO 10 M= 1 , MM 

D I ST=QUAD (M ) * ( DSQRT ( A+B ) -DSQRT ( A - B ) )* . 5 
. + (DSQRT ( A+B ) +DSQRT ( A - B ) ) * . 5 
THCOS= ( A -DI ST**2 ) / B 
10 RIGHTl=RIGHTl+6./DFL0AT(MM)*RH0 

. /DSQRT ( ( DSQRT (A-B)+DIST) * ( DSQRT ( B+A ) +DI ST) ) 

. * (R J -RHO*THCOS ) /DI ST 

. * ( Z J - ZETA ( 2 ) -DZ 1DR 1* (R J*THCOS -RHO) ) /DI ST 
.-•'•-(UR 

. * ( R J*THCOS -RHO ) + ( Z J - ZETA ( 2 ) ) *UZ ) 

RETURN 

END 

FUNCTION RIGHT2(RH0) 

IMPLICIT REAL* 8 (A-H.O-Z) 

REAL*8 QUAD (30) ,VECRH0(3) ,ZETA(3) ,C(30,3) ,R(30) ,ORD(30) 

LOGICAL UPFLOW 
COMMON 

./B/C.R.ORD.BPAR 
. /D/SMALL, P 

. /C/QUAD .RJ.ZJ.R2. MM . NR . UPFLOW 
EXTERNAL EXPIN 

IF (RHO.GT. (R( 1)+SMALL) ) GOTO 998 
VECRH0(1)=R(1) 

VECRH0(2 )=R ( 1 )+SMALL 
VECRHO ( 3 )=R ( 1 )+SMALL*2 

CALL ICSEVU(R,0RD,NR+3,C, 30, VECRHO, ZETA, 3, IER) 

DER 1= (ZETA ( 2 ) - ZETA ( 1 ) ) / SMALL 
DER2= ( ZETA ( 3 )+ZETA ( 1 ) -2*ZETA (2 ) ) /SMALL**2 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C ZETA(2) WILL NOW REPRESENT ZETA IN. [ZER0PL,R(1) ] 

C DZ1DR1 WILL NOW REPRESENT DERIV IN [ZER0PL,R(1) ] 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

ZETA ( 2 )=ZETA ( 1 )+ (RHO-R ( 1 ) )*DER 1+ (RHO-R ( 1 ) )**2/2 . *DER2 
DZ1DR1=DER1+(RH0-R( 1 ) )*DER2 
GO TO 997 

998 VECRHO (1)=RH0-SMALL 
VECRHO ( 2 )=RHO 
VECRHO ( 3 )=RHO+SM ALL 

CALL ICSEVU (R , ORD , NR+3 , C , 30 , VECRHO , ZETA , 3 , IER) 

DZ1DR1=(ZETA(3) -ZETA( 1 ) )/2 . /SMALL 
997 W= ( ZETA ( 2 ) +DSQRT ( ZETA ( 2 ) ** 2+RHO** 2 ) ) / P 
S= ( - ZETA ( 2 ) +DSQRT ( ZETA ( 2 ) ** 2+RHO**2 ) ) / P 
IF (S.LE.O) S=0. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C FLOW DOWN 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

IF (UPFLOW) GO TO 811 


UZ=- ( - 1 . +EXPIN (R2*V ) /EXPIN ( R2 )+ 

. (DEXP( -R2) -DEXP( -R2-W) )/EXPIN'(R2)/R2/ (W+S) ) 

UR=- ( ( 1 . -EXPIN (R2*W) /EXPIN (R2 ) )*DSQRT(S/ (S+W) ) ) 

GO TO 812 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C FLOW UP 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

811 UZ=- (1. -EXPIN (R2*W) /EXPIN (R2)- 

. (DEXP(-R2) -DEXP(-R2*W) ) /EXPIN (R2)/R2/ (W+S) ) 

UR=- ( - ( 1 . -EXPIN (R2---W) /EXPIN (R2 ) )*DSQRT(S/ (S+W) ) ) 

812 B=2*RHO*RJ 

A=RHO** 2+R J**2+ ( ZETA ( 2 ) - Z J )' **2 
RIGHT2=0 . 

DO 10 M= 1 , MM 

DIST=QUAD(M)*(DSQRT(A+B)-DSQRT(A-B))*.5 
. + (DSQRT (A+B I+DSQRT (A-B ) )* . 5 
THCOS= ( A -D I ST**2 ) / B 
10 RIGHT2=RIGHT2+6./DFL0AT(MM)*RH0 

. /DSQRT ( (DSQRT(A-B)+DIST)- V (DSQRT(B+A)+DIST) ) 

.*(ZJ-ZETA(2))/DIST 

.*(ZJ-ZETA(2)-DZ1DR1*(RJ' V THC0S-RH0))/DIST 

.*(UR 

. * (RJ*THCOS-RHO)+ ( Z J - ZETA ( 2 ) ) *UZ ) 

RETURN 

END 
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Abstract 

In order to determine a protein's tertiary structure, large, well ordered single crystals are 
required for x-ray analysis. Producing such crystals is often the rate-limiting step because many 
protein crystals grow slowly and they often reach a terminal size which is too small to be useful for 
x-ray diffraction studies. In this paper, we present our study of several mechanisms which may 
reduce crystal growth rates and/or terminate crystal growth entirely. On the basis of our analysis, 
we find that salt gradients which change the local chemical potential of the protein are insufficient to 
account for the slow crystal growth rates which have been reported. Contaminants which adsorb 
protein from solution may reduce the effective protein concentration, but the impurity’s 
concentration and its affinity for protein are unknown. Association of protein molecules in bulk 
solution can reduce the monomer concentration significantly, but extant theory and experiment are 
not sensitive enough to determine the actual concentration of aggregates in solution. For systems 
of interest, shear-induced effects were found to be too weak to interfere with normal binding of 
incoming protein molecules. Although we found that most crystal growth occurs in a regime where 
both interfacia] kinetics and diffusion influence crystal growth, the role of mass transfer rates on the 
terminal size of crystals is unknown, primarily because no data exist which cover the size range of 
interest (0.1mm - 1mm in length). 

Experimental studies of growth of large protein crystals are essential if mechanisms by which 

crystals stop growing are to be elucidated. Growth rate measurements for a wide range of crystal 

\ 

sizes, coupled with measurements of system properties which vary with time, may reveal the 
factors responsible for this puzzling behavior. Several hypotheses have been advanced to explain 
growth cessation, but none have been verified by experiment. For example, if termination of 


growth occurs by accumulation of defects, x-ray studies of the crystals at different stages in their 
growth should reveal some qualitative changes in the nature of the protein crystal which is related to 
its growth behavior. 

Introduction 

Protein crystals are useful for determining the tertiary structure of biological molecules but the 
task of preparing crystals is difficult because many proteins do not form crystals readily or form 
crystals which are unsuitable for x-ray analysis. Protein crystals differ from more f amili ar 
inorganic crystals in several ways. In many systems, crystals tend to grow at a constant rate, up to 
a point, and then stop growing. The molecules are held in the crystal lattice by weak hydrogen 
bonds (AH = -3 kcal/mole to -6 kcal/mole in vacuo ) which are not well characterized [Creighton 
1983] and each molecule has relatively few contacts with its neighbors [McPherson 1982], 
Additionally, the bonding structure and the physical nature of the crystal are not well understood. 
Moreover, the cry'tallization process itself is complicated by small changes in temperature, pH, or 
the precipitating agent concentration. Although it has been shown that the precipitant concentration 
inside the crystal equilibrates with the external solvent [Wyckoff 1986], the concentration of 
precipitants in crystals has not been measured. Thus the mechanism by which precipitants act and 
their role inside the crystal are unsolved mysteries. 

Some of the difficulty of understanding protein crystallization arises from the complexity of a 
system containing protein, precipitating agent, buffer, and solvent Interactions between these 
components makes it nearly impossible to predict whether a given set of conditions is suitable for 
protein crystal growth. Thus, guidelines for producing protein crystals are largely rules of thumb. 
Results from groups working on the same model system under similar conditions can differ 
significantly, as evidenced by disparities in the solubility of lysozyme reported at 20*C in 
5%(wA0 NaCl at pH 4. Values range from 1.7 mg/ml [Pusey et al. 1986] to 4.3 mg/ml [Kam et 
al. 1978, Feher and Kam 1985] with an intermediate value of 3.5 mg/ml [Fiddis et al. 1978]. 


When used to determine the concentration dependence of the crystal growth rate these solubilities 
lead to different kinetic expressions which makes it difficult to decide on the controlling 
mechanism(s). 

As though studying crystal growth by diffusion of protein from the bulk were not complicated 
enough, fluid flow may also influence matters. It is known that forced convection in the form of 
stirring tends to produce large numbers of small crystals rather than the large single crystals 
desired. Recently it was suggested that natural convection engendered by the protein concentration 
gradient near the crystal surface disrupts the orderly deposition of protein molecules [Bugg et al. 
1984]. This fluid motion might: (i) promote rapid transfer where incoming protein molecules are 
forced to bind at a random site rather than the specific sites that lead to well ordered crystals, or 
(ii) alter the state of the protein at or near the surface surface through fluid shear. In the first case 
the convective mass transfer rate is envisioned as being faster than the attachment rate so defects 
form due to improperly bound or misoriented molecules. Accordingly, growth ceases due to an 
accumulation of ?.?rors when the surface concentration of defects reaches the point where additional 
molecules cannot find a suitable attachment site [Kam et al. 1978, Feher and Kam 1985] . This 
implies that neighboring molecules on the surface are misaligned and the degree of local disorder is 
extremely high or that binding requires the cooperation of a large number of properly oriented 
molecules on the surface. Unless these defects arise suddenly, crystal disorder would increase 
gradually moving towards the crystal surface (with a corresponding decrease in the x-ray 
resolution). According to the second mechanism, shear stresses at the crystal surface strip the 
molecules from the crystal surface or align the protein molecules so that molecules near the surface 
cannot assume the proper orientation for binding to the crystal. The shear stress on the protein 
molecule would need to be strong enough to denature the protein molecule by breaking its internal 
hydrogen bonds and unfolding the molecule. 

Finally, contaminants in the solution may accumulate at the surface of the crystal and either 
poison the surface so that no further addition is possible, or induce some change in the protein near 


the crystal surface which prevents it from binding. 

In this paper, two main questions are addressed: (i) Why do protein crystals grow so slowly? 
and (ii) Why do protein crystals reach a terminal size? To answer the first question, we have 
examined three mechanisms which might limit crystal growth rates: salt rejection at the crystal-fluid 
interface, contamination by impurities which adsorb protein, and the aggregation of protein 
monomer in the bulk (which would lower the effective concentration of protein available for crystal 
growth). Our analysis indicates that any salt gradients due to salt rejection are too small to produce 
the effects needed to reduce the diffusion-limited crystal growth rates to measured levels. The 
presence of dust or other contaminants may significantly reduce the effective protein concentration 
if the protein adsorbs strongly to the dust surface. Parameters which affect the importance of 
surface adsorption, such as the protein's affinity for the dust and the dust concentration found in 
typical protein solutions, are not currently knowri. Aggregation of protein may significantly reduce 
the amount of protein monomer in solution, but both the measurements and the current theory for 
determining size distributions are inadequate to determine the actual state of aggiegation of the 
protein in solution. Further revisions of the theory and a more comprehensive set of experiments 
may provide enough information to determine whether protein aggregation plays a role in limiting 
crystal growth rates. 

Several hypotheses concerning the effect of fluid flow on the terminal size of protein crystals 
were investigated by examining several mechanisms in which fluid shear at the crystal interface 
interferes with normal addition of protein molecules to the surface and by performing a 
quasi-steady analysis of the mass transfer to the crystal surface. The effects of shear due to fluid 
flow were found to be several orders of magnitude too weak to disrupt proper protein attachment to 
the crystal surface or to align the protein molecules in the vicinity of the surface. An analysis of 
mass transfer indicates that crystal growth occurs in a regime where diffusion and the attachment 
reaction both play roles in limiting the crystal growth rate. 

Formulation of a a theory to explain the peculiarities of protein crystal growth is hindered by 


the lack of reliable data on growth rates of large crystals. Extant results were obtained with crystals 
no larger than approximately 100 mm in length and size effects are not yet apparent in these 
crystals, so vital information about the slowing growth rate is missing. A thorough study should 
be performed wherein larger crystals are produced in order to relate growth rates to system 
properties which may change over the course of the crystal growth process. 

Salt Rejection and Protein Crystal Growth 

There are at least two possible ways in which the rejection of salt (precipitating agent) at the 
interface may influence protein crystal growth: (i) a "blowing" velocity away from the crystal 
surface which would slow the transport of protein to the crystal; and (ii) the alteration of the local 
protein solubility which would reduce the driving force for diffusion. In the first case, the blowing 
would appear in the "crystallization" flow which arises from the diffusion of protein to the crystal 
surface. This crystallization flow is related to the growth rate of the crystal, dR/dt, by 

v^n-S^O-p^pp (1) 

where v f is the fluid velocity at the interface, n is the unit normal directed outward from the crystal 
surface, p c is the crystal mass density and p f is the fluid density at the interface. The crystallization 
flow is directed towards the crystal surface if p c > p f , and away from the crystal if p c < p f . A rising 
convective plume has been observed [Pusey 1986, personal conversation], indicating that the fluid 
at the interface is less dense than the bulk fluid. Because the crystal is denser than the bulk fluid, 
p c > p { and the rejection of salt is insufficient to induce a blowing away from the crystal surface. 
Earlier calculations [Grant 1985] showed that the convective protein flux due to crystallization flow 
is approximately 1 % of the diffusive flux for a spherical protein crystal. The effect of salt rejection 
is to reduce this small convective flux, and can be neglected without appreciable error. 

The effect of variations in the local protein solubility can be examinied by considering the 
growth of a spherical protein crystal under diffusion control. At the crystal surface, salt is rejected 


and must diffuse to the bulk solution so that a salt concentration gradient exists along with the 
resulting gradient in the local solubility of the protein. In the diffusion limit, protein in the liquid at 
the ciystal interface is at the solubility concentration and is in equilibrium with the protein in the 
crystal. Since the crystal form is reported to be insensitive to the salt concentration, it may be 
reasonable to assume the chemical potential of the crystalline protein is independent of the salt 
concentration. In this case, the chemical potentials of the protein in solution and in the crystal are 
equal wherever the protein concentration equals the local protein solubility because we could, in 
principle, place a crystal in a saturated protein solution without producing a change in either the 
crystal size or the protein concentration of the solution. The chemical potential of the protein 
depends on the other species which are in solution, so the protein's chemical potential is no longer 
directly proportional to its concentration. The salt concentration gradients due to rejection at the 
crystal interface may alter the gradients in the protein's chemical potential with the result that the 
flux of protein to the crystal surface is less than the flux one would expect from examining only the 
protein concentration gradient. 

To set this out in mathematical form, we first express the flux of protein (species 1) in terms of 
its chemical potential [Cussler 1982] 

where (ij is the chemical potential of the protein, D 0 = thermodynamic diffusion coefficient, and 
Cj= molar concentration of protein. The chemical potential of the protein is given by 

p^po + kTlnx^ (3) 

where Pj° = standard state chemical potential, Xj= mole fraction of protein, and Yj=activity 
coefficient of the protein. Since the chemical potential of the protein at the solubility concentration 
is constant, it follows that |ij so1 =p 1 °+ kT In x^ 1 Yj so1 = constant. The flux relation given by 
Equation (2) is unchanged by adding the gradient of a constant, so the flux can also be expressed as 


(4) 


j -i^ C > V( VO=-T^ c . V1 " 
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One form for the activity coefficient is [Debenedetti 1984] 


-Kx, 




(5) 


where is the activity coefficient of protein when in the limit of an infinitely dilute solution. If K 
is a constant independent of salt concentration and varies with salt concentration, then the 
chemical potential of the dissolved protein can always be written in terms of the local Yj“: 


|i° + kT In Xj Yj = |i° + kT In 
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( 6 ) 


Although this manner of adjusting Yj“ to satisfy the solubility constraints is ad hoc, it makes the 
mathematics somewhat simpler by absorbing the salt dependence into Xj so1 so that substitution of 
equations (5) and ( 6 ) into Equation (4) gives 


J = - 


D 0 Cj V In 


V, *0l v 

-K(Xj - Xj ) 


sol 


Li 


(7) 


A quasi-steady mass balance around a growing spherical crystal is 

( R f-) c x= r2D o c ii to 


X. -K(x | ■ Xj ) 

"soT e 
X 1 


( 8 ) 


where C x is the molar concentration of protein in the crystal. If the following dimensionless 
variables are introduced: 
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( 9 ) 


the mass balance becomes 


x >di ln 


X = R/R 0 

z = R/ r = R Q X./r 

^•Do/rJ 


V! 501 X 

Xj -K(Xj - Xj ) 

sol 

X 1 


= -Pe (C x /Cj) 


Pe = X-^7 = -&-dE. 
dx D 0 dt 


( 10 ) 


where the relation Cj = Cj-Xj has been used and Cp will be taken as 55 molar. Transforming the 
equation from an expression for the activity into a differential equation for Xj yields: 

sol 

sol Q2 

rr tc \ -u r^v _v / v ^ 

dx. 


k l _ 


Pe (Cjj/Cr) + (Kxj-x^x” 1 )^- 


dz 


1 -Kx 


( 11 ) 


The salt concentration was calculated from the quasi-steady diffusion profile 

Csalt = C sa i t tC » + ^Qalt z 


( 12 ) 


where C sa]lo = 50 mg/ml, and AC S = C saltsurf - C sa]Uoo . The solubility was assumed to obey the 
expression [Feher and Kam 1985] 


sol 

Xj =ae 


salt 


Combining equations (12) and (13) gives Xj so1 as a function of position: 


sol sol -B ACeoUZ 
X 1 =x l,~ e 


(13) 


( 14 ) 


Equation (14) was substituted into (1 1) to yield 
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( 15 ) 


I 


, -PAC , z 

PeCC^^-fPAC^x^PAC^Kx^^e sal 
dz 1 - KXj 


and the initial value problem could then be integrated numerically from the crystal surface (z=l , 
x ,=xj s so1 ) to infinity (z=0, x } =Xj < J. The resulting value of Xj ^ can in turn be used to calculate 
the nominal growth rate which would be predicted from the diffusion limit in a uniform salt 
concentration field, viz. 


Pe 


nom 


= £l 
c x 



sol v 

' x i ) 
1 


(16) 


and the ratio of the actual and nominal growth rates was determined as a function of system 
properties. 

For a dilute protein solution, Kxj«l [Debenedetti 1984], and can be neglected in the 
denominator of Equation (15), while the Kxj term in the numerator will be neglected in order to 
obtain an upper bound on the effect of salt rejection on he diffusion rate. In this case. Equation 
(15) can be integrated analytically to yield 


sol 

x = x e 

1 1 oo 
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(17) 


The value of Xj ^ obtained from Equation (17) can be substituted into (16) to obtain 


Pe = 


pAC 


salt 


^PAC S alt _ j 


Pe 

nom 


(18) 


Although the form of the relation between the nominal Peclet number and the actual Peclet 
number has been established, the apparent reduction in diffusion rate depends on AC salt , which has 
not yet been determined. Recall that the quasi-steady salt concentration profile, Equation (12), 
was used to obtain Equation (18), but the surface concentration was left unspecified. A mass 
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balance on the salt rejected at the crystal surface yields: 




= D«, 


dQalt 


'salt,sur f ^salt “3F 


(19) 


where s = segregation coefficient of the salt (C salucrysla] / C sall surface ); 0 < s < 1. In terms of the 
dimensionless variables introduced in Equation (9), Equation (19) is 


Pe (D Q / E^ all ) (1 s) C saW - 


( 20 ) 


so that 


AC = C — C = C 

salt salt,surf ^sali,®® salt,' 


1 


l-Pe(D n /D salt )(l-s) 


( 21 ) 


For proteins, D 0 / D sa]t «* 10 1 , and typical Peclet numbers for lysozyme crystals are less than 6 
x 10'^. The maximum concentration difference occurs when s = 0 (total segregation of the salt), in 
which case AC sa]t < 6 x 10' 3 C sa]Uoo . At a bulk concentration of 50 mg/ml NaCl, this gives 
AC sa] ,< 0.3 mg/ml. For the actual growth rate to equal one tenth the nominal growth rate, p should 
be approximately 10, but reported values of lysozyme solubility put the value of P closer to 10 1 
[Ataka and Tanaka 1986]. Thus, Pe = 0.985 Pe nom , and salt rejection cannot reduce the diffusion 
rate to the point where crystal growth is entirely diffusion controlled. 

If the protein molecules are treated as large ions and the protein activity coefficient obeys the 
Debye-Hiickel law, the chemical potential of the protein molecules decreases as the ionic strength 
increases. This leads to a greater driving force for diffusion and a corresponding increase in 
growth rate, a trend which would make diffusion less likely to limit the crystal growth process. 

The solutions from which protein crystals are usually grown, however, are too concentrated in salt 
for the limiting form of the Debye-Hiickel law to hold and it is unreasonable to expect the protein 
molecules to behave like simple ions [Tanford 1961]. The actual form of the activity coefficient, 
then, may partially acount for the slow growth of protein crystals. Any mechanism which requires 


large gradients in salt concentration due to salt rejection cannot account for the slow growth rate 
unless the protein is extremely sensitive to small changes in salt concentration. 


Effect of Impurities on Protein Crystal Growth 

The presence of dust or other impurities affects protein crystal growth in two different ways: 1) 
the dust surface acts as a nucleation site for the protein; and 2) the protein molecules may adsorb to 
the surface of the dust particle, thereby reducing the effective protein concentration in the solution. 
The nature of the dust surface is too ill-defined to determine quantitatively how the nucleation rate is 
affected, so we will address the effect of protein adsorption in this section. 

Suppose that there are a certain number of sites, n T , on the dust surface to which protein 
molecules will bind. Further, assume that the rate at which the protein adsorbs is proportional to 
the product of the protein concentration and the concentration of available sites, and that the 
desorption rate is proportional to the concentration of occupied sites: 

TT=- c J k «i c p n *- k i <n T- n ‘ )) ( 22 > 


where C p =molar concentration of protein, C d =molar concentration of dust, n*= number of 
available binding sites/dust particle, k ad =adsorption rate constant, and k d =desorption rate constant. 
At equilibrium, dn*/dt=0, and we have 


_ k ad _ ftp - n* 

k d n * c „ 
d p 


0eq 

(1 ~0eq) Cp >e q 


(23) 


where 0 is the fractional coverage of the particle surface. 

A mass balance on the protein in solution yields: 

C p = C P,0-"T C d e < 24 > 

where C p 0 is the original protein concentration in solution. Equation (23) can be substituted into 


(24) so that 


C D = C nn -n_ C 


KCp^ 


p .“l M T-d , +KC 


P,eq 


(25) 


which can be solved to yield: 


C P.eq — 


- (1 - KC p 0 + Hj. C^) + Jq - KC p 0 + Or C^) 2 + 4 KC p 0 


(26) 


2K 


We can take a typical dust particle to be approximately 0.5 pm in diameter, and if the particle is • 

assumed to be completely covered by a monolayer of protein molecules (taken to be cubes with a 
length of twice the hydrodynamic radius of 20 A), approximately 5 x 10 4 protein molecules can be 
adsorbed on each particle of dust. This indicates that a 5%(w/v) solution of lysozyme can be • 

completely depleted by adsorption if ^>6.9 x 10' 8 mol/liter, provided that K»l. 

Contamination by dust, then, can serve to reduce the effective protein concentration by a large 
amount. As the protein is depleted by the growing crystals, the molecules on the dust can desorb • 

and replenish the protein concentration in solution. This would maintain the protein level in the 
solution at a fairly constant level, which might explain the nearly constant rates of growth which 
have been observed. • 

Aggregation of Protein 

The state of aggregation of protein in solution affects the crystal growth rate by changing the • 

rate at which protein diffuses to the crystal surface and reduces the effective concentration of 
protein available for growth. The diffusion coefficient of a particle is inversely proportional to its 
hydrodynamic radius, so that the aggregates have smaller diffusion coefficients than monomers and • 

will diffuse to the surface slower. Additionally, if crystal growth occurs by addition of aggregates 
of specific sizes (e.g. monomers, dimers, or trimers), then the formation of non-participating 
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aggregates acts as a protein sink, thereby reducing the growth rate. Kam et al. [1978] and Feher 
and Kam [1985] have proposed a method for determining a particle size distribution based on 
quasi-elastic light scattering measurements and a thermodynamic model of nucleation. The 
discussion below is in two parts: (i) an analysis of the experimental measurements, and (ii) a 
description of their model. 

According to Kam el al. [1978] and Feher and Kam [1985], a solution illuminated with 
incident light of frequency v o will scatter light with a power spectrum given by: 



A Vj 


( V s" V c/ + < AV p 2 


(27) 


where j = number of molecules in an aggregate, v s = frequency of scattered light, AVj = 

(2Djq 2 /2n), q = scattering vector, Dj = diffusion coefficient of a j-mer, and Cj is proportional to the 
number density of the j-mers. The power spectrum obtained experimentally is usually 
approximated by a single Lorentzian: 


S(v) 


A 

(v - v o ) 2 + (Av) 2 


(28) 


where A is a constant and Av is the effective linewidth of the power spectrum. The parameters A 
and Av can be found from a non-linear least-squares fit of the spectrum to the form given by (28) 
or by transforming the data into the form 

Y(v) = 1/S(v) = (V s ~ - V fl )2 = mx + b (29) 

A A 


where m = A' 1 , x = (v $ - v o ) 2 , and b = (Av) 2 /A. In these coordinates, Av = (b/m) 1/2 and the 
parameters b and m can be obtained from a linear least-squares fit of Y(v). The linewidth is a 
rough measure of the state of aggregation of the protein because it has a maximum value when the 
protein is entirely monomeric, and decreases as the fraction of aggregates in solution increases. 


The amplitude, A, also contains information on the concentration of aggregates: as the number of 
aggregates increases, so does the amplitude of the power spectrum. Analysis of experimental data 
in this manner yields only an "average" degree of aggregation rather than information on the actual 
size distribution of the protein. Although some sort of average diffusion coefficient can be 
determined direcdy from the linewidth, more detailed information on the particle size distribution 
cannot be obtained from the scattering spectrum alone. For this reason, Kam el al. chose to 
measure how the normalized linewidth, Av*=Av/AVj, changes with protein concentration and to fit 
the parameters of their model so that the agreement between theory and experiment is maximized. 
The resulting particle size distribution, which is consistent with both theory and experiment, may 
not accurately represent the true size distribution in solution. As we will show later, the proposed 
model is not sensitive enough to discern the true distribution unambiguously. 

The model is based on the thermodynamics of aggregation of monomers to form nuclei of 
different sizes. The change in free energy required to form a nucleus from monomers is usually 
considered to consist of two contributions: a negative bulk contribution which acts to promote 
aggregation and a positive surface component which arises from the need to support the additional 
surface area of the aggregate. The formation of small nuclei is thermodynamically unfavorable 
because the free energy change is positive due to the large surface to volume ratio of the aggregates. 
As the nucleus grows, the addition of monomers becomes less unfavorable until the nucleus 
reaches a critical size beyond which further addition of monomers is favorable. Any nucleus which 
exceeds the critical size will grow spontaneously until it reaches macroscopic size. The critical size 
of the nucleus depends on how the surface area and volume of the nucleus varies with aggregate 
size. For a spherical nucleus, the volume (bulk) is proportional to the number of molecules, j, in 
the aggregate, while the surface area is proportional to j 2 ^ 3 . The standard free energy change for 
the formation of an aggregate from j monomers is 

AGj^O-DGb + PO^-DGs (30) 

where G B = bulk free energy per molecule (assumed constant), G s = surface free energy per 

- 14 - 


molecule (also constant), and (3 = geometrical factor which gives the surface area per molecule. 
The standard state conditions are 1 %(w/v) protein, at the temperature, pH, and salt concentration 
of the system. The actual change in Gibbs free energy is then 


AGj = (j - 1) G fi + P (j 2/3 - 1) G s - (j - 1) kT In Cj, (31) 

where Cj. is the total protein concentration expressed in units of %(w/v). Differentiation of (31) 
shows that AGj reaches a maximum at the critical size 

1 " fPGs/lcT (3 2 ) 
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The equilibrium constants for the aggregation reactions 
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(34) 


Although this approximation is reasonable for large clusters, it becomes less accurate for s mall 
aggregates such as those which may be present during the pre-nucleation stage. If the expression 
for AGj 0 from Equation (30) is substituted into Equation (34), the result is 


In K. = - 

j 



(35) 


Evaluating Equation (35) for j=l and j— »«> shows that the equilibrium constants are related to die 
free energy terms as follows: 


(36) 
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Substitution of (36) into (32) leads to 
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and the free energy barrier for nucleation is 
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Kam et al. chose to use Kj and K eo /K 1 as the adjustable parameters in their model. 

When AGj*/kT » 1, the formation of critical-sized nuclei is slow, and the system can reach a 
state of "quasi-equilibrium” in which particles smaller than j* are in equilibrium with each other and 
the concentration of larger particles is negligible. The concentration of aggregates during the 
pre-nucleation stage was calculated from the quasi-equilibrium approximation: 


C. =K . C. . C. 
j j-i j-i i 

C. =0 

J 

2jS =c t 


; 
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In order to calculate the quasi-equilibrium distribution, all that is required is to select different 
values of Kj and K^, calculate j* from Equation (37), and solve Equation (39) for the aggregate 
concentrations. The distributions obtained from the reported values of Kj= 0.065 %( w/v)' 1 and 
K^/Kjs 35 are shown in Figure 1 for different protein concentrations. 

Once the quasi-equilibrium distribution is known, the approximate power spectrum can be 
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calculated by substituting the Cj into Equation (27) and truncating the summation at j*. Kam et al. 
then plotted the normalized linewidth, Av*, against the normalized protein concentration, KjCp for 
different values of K oo /K 1 (solid curves in Figure 2). The normalized linewidths obtained from 
experiment were plotted against protein concentration on the same semilogarithmic scale, and then 
shifted until agreement between the experimental data and one of the theoretical curves was 
attained. The value of Kj can be calculated from the shift required to align the experimental and 
theoretical curves, while the appropriate value of K oe /K 1 is that corresponding to the theoretical 
curve. Unfortunately, imposing this model on the system can lead to erroneous interpretation of 
the data. To illustrate the problem, we have created a sample particle size distribution in which the 
aggregates cannot grow larger than j=4 and the equilibrium constants are Kj= 1 %(wA')* 1 , 

1^= 5%(w/v)' 1 , and K 3 = 3%(w/v)' 1 . The power spectrum was calculated as a function of 
concentration according to Equation (27) and the linewidth was calculated from Equation (29). We 
then obtained an approximate fit according to Kam's method described above (Figure 2) and found 
that the best-fit parameters are Kj= l%(w/v) -1 and K oe /K 1 = 130. The disparity between the trial 
distribution and the best-fit distribution is shown by comparing the fraction of monomer at 
l%(w/v) protein concentration: the trial distribution contains approximately 26% monomer while 
the best-fit distribution contains essentially no monomer. This method of determining the size 
distribution is simply not sensitive enough to differentiate among all the possible distributions 
This insensitivity to the actual particle size distribution is not the only problem with Kam's 
model; the quasi-equilibrium approximation breaks down at concentrations in the range of interest 
to crystallographers. The quasi-equilibrium approximation is valid only when the free energy 
barrier is large compared with kT. At higher concentrations, the free energy barrier is relatively 
low (Figure 3) and critical nuclei form rapidly enough that their concentration cannot be neglected. 
When the quasi-equilibrium approximation is applied in these circumstances, the protein is forced 
to distribute itself among an artificially small number of aggregate sizes with the result that the 
model predicts large "jumps" in the mass fraction of monomer due to small changes in total protein 


concentration. At concentrations where small changes in protein concentration decrease j* by one, 
the fraction of monomer increases discontinuously (Figure 4a). The oscillations seen in the 
normalized linewidth (Figure 4b) are a consequence of these monomer jumps. The curves in 
Figure 2 do not show this behavior because they were terminated at concentrations below those at 
which Av* starts to increase. 

The thermodynamic model of aggregation presented here is based on two assumptions that can 
have large effects on the predicted particel size distribution. The first is that the protein solutions 
behave as ideal solutions. Even though the solutions are dilute, they are highly supersaturated and 
could deviate significantly from ideal solution behavior. The second assumption is that both the 
bulk free energy and the surface free energy remain constant with crystal size. At least in the case 
of the surface energy, it would be reasonable to expect some size dependence because it makes little 
sense to speak of bulk and surface properties when the aggregates are as small as four or five 
molecules. The most likely possibility is that the surface energy is small for small aggregates and 
gradually increases until it reaches its macroscopic value. This would make the formation of small 
aggregates more favorable and reduce the monomer concentration. 

On the basis of the analysis of Kam's work, aggregation of protein cannot yet be eliminated as 
a mechanism which slows protein crystal growth. There is enough ambiguity in the experimental 
results that the amount of protein in aggegates cannot be determined with any confidence. 
Definitive studies of particle size distributions or, at the least, average particle sizes under growth 
conditions may show whether or not protein aggregation plays an important role in controlling 
crystal growth. In addition, a refinement of the theory to include size dependent effects may show 
that the effective concentration of monomer is much lower than curently believed. 

Effects of Fluid Flow on Protein Crystal Growth 

Bugg et al. [1984] suggested that the terminal size of the crystals could result from the size 
dependence of the natural convection on crystal size. If the force required to break the crystal 


bonds is comparable to that of the shear stress at the crystal surface due to natural convection, then 
the flow field may be strong enough to strip protein molecules from the surface. Because the shear 
stress at the surface increases with crystal size, this would act as a self-limiting process. We can 
determine if this is a likely explanation for the phenomenon by performing an order of magnitude 
estimate of the viscous stress based on the velocity profile for a vertical flat plate. 

From the velocity profile presented by Schlichting [1979], the shear stress is given by 



where: t=shear stress, p=fluid viscosity, r=shearrate, g=gravitational acceleration, Ap=density 
difference, p M =bulk fluid density, v=kinematic viscosity, R=crystal radius. According to Equation 
(1.1), a 1mm diameter spherical crystal in a 5%(w/v) lysozyme solution under unit gravity would 
experience a shear rate on the order of 100 s' 1 . This analysis ignores the effect of the Schmidt 
number, Sc=v/D, on reducing the shear rate. At high values of Sc, such as Sc=10 4 for lysozyme, 
T would be much smaller than the value obtained above [Schlichting 1979]. Nevertheless, if we 
use r=100 s* 1 , any effect of fluid flow on protein crystal growth will not be prematurely ruled out. 
The shear stress acting on the crystal surface in a solution where p = 1 x 10* 3 Pa-s is approximately 
0.1 Pa. If we take Fiddis* [1978] approximation of the lysozyme molecule being a cube 30.9 A on 
a side, then the shear force acting on the molecule is approximately 10' 18 N. 

We can take Fiddis’ [1978] value of the surface energy (7.5 kJ/mol) as representative of the 
strength of the crystal bond, but the form of the potential is still unknown. For want of a better 
estimate, we approximate the potential as a Lennard- Jones 6-12 potential with electrostatic 
interaction: 


where A and B are constants for the particular bond of interest, and q 2 are the partial charges on 
the atoms in the bond, e 0 is the permittivity of free space, and r is the separation between the atoms. 
Hagler et al. [1974] determined the values of A, B, qjand q 2 for the bonds in amide crystals and 
found that the form of Equation (1.2) accounted adequately for the observed interactions. From the 
values of A, B, and q reported by Hagler et al., the bond energies were calculated from the value of 
U at the equilibrium separation (where F = -dU/dr = 0) and the maximum attractive force was 
calculated by determining the force where dF/dr = -d 2 U/dr 2 = 0. The bond energies ranged from 
-0.3 kJ/mole for the bond between the carbonyl carbon and the amino hydrogen, to -29.5 kJ/mole 
for the bond between the carbonyl carbon and carbonyl oxygen. The bond with an energy closest 
to the 7.5 kJ/mole reported by Fiddis is that between hydrogen and nitrogen, which has a strength 
of 5.6 kJ/mole and requires a maximum force of 2.2 x 10' 11 N per bond to break. If the 
electrostatic potential is neglected, then the bond with the lowest breaking force is that between a 
non-carbonyl carbon and the amino hydrogen, which requires F=8.3 x 10 13 N/bond to break. 

The force generated by free convection is approximately six orders of magnitude too small to strip 
molecules from the crystal surface. 

Even though the shear stess cannot remove molecules from the surface, it may impart some 
preferred orientation to the molecules near the surface so that they are unable to Fmd the proper 
alignment for addition to the crystal. To test this hypothesis, we compare the characteristic rates of 
the processes: if the velocity gradient at the surface tends to align the protein molecules in a 
preferred orientation, it acts at a rate which is comparable to the shear rate, T; the characteristic 
rotational rate is given by the rotational diffusion coefficient, D rot = kT/SttpR 3 , where R is the 
hydrodynamic radius of the protein [Tanford 1961]. The hydrodynamic radius of lysozyme is 
approximately 20 A [Tanford 1961], so D rot = 2 x 10 7 s' 1 while the shear rate is certainly less than 
100 s* 1 , so the ratio of rotational and shear rates is 2 x 10 5 . This indicates that randomization of 
the protein molecules occurs much faster than any orientation imposed by the shear flow. Once the 
molecule has reached the surface and formed some sort of bond, however, the situation is different 


because the molecule is constrained at the binding site. If an improperly oriented protein molecule 
bonded to the surface, it would be possible for other protein molecules to orient themselves 
properly with respect to the first molecule. In this manner, the crystal would be made up of regions 
with the same local orientation and the average ordering for the crystal as a whole would be 
reduced. If the crystal does not consist of this ensemble of regions, it would indicate that some 
ordering takes place, possibly due to the electrostatic interactions between charged groups. 

The effects of electrostatic interactions on the interatomic potential were mentioned only briefly 
above, but they exhibit a strong influence on the nature of the bonds. The electrostatic contribution 
is large compared with the other components of the hydrogen bond, so it is important to know the 
length scale over which it acts. The natural length scale for electrostatics is the Debye shielding 
length, K" : 
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where e=dielectric constant of water, k=Boltzmann's constant, T= absolute temperature, 
N A =Avogadro's number, Z=charge on ion, e=electron charge, and I=ionic strength. The Debye 
length is approximately 3 A for 50 mg/ml NaCl, which indicates that the effects of the electrostatic 
potential are substantially reduced when the atoms are separated by distances greater than the typical 
bond lengths reported by Hagler et al. Any orientation due to electrostatic interactions, therefore, 
would occur when the incoming protein molecule is practically bonded to the crystal surface. 

The denaturation of proteins by the shear field was suggested by Marc Pusey in a conversation 
with M. L. Grant. The following model, suggested by W. B. Russel of Princeton University, was 
used to investigate the possibility of shear-induced denaturation. If the protein molecule must be in 
a particular conformation in order to bind to the crystal surface, crystal growth may be hindered if 
the shear stress due to fluid flow is sufficient to change the protein's conformation. For this 
analysis, consider a molecule of protein to be spherical as shown in Figure 5. The molecule is 


maintained in this conformation by a single hydrogen bond placed at point A and is "hinged" at 
point B. If the shear forces on the molecule are sufficient to break the bond at A, the molecule will 
open up in the yz plane and the molecule will be unable to bond to the surface. The force on the 
bond at A can be determined by calculating the torque about B due to creeping flow past the sphere 
and determining the equivalent force to place at A. 

Over any element area on the molecular surface, the magnitude of the torque is given by 
but only the component in the ±x direction contributes to opening the hinge in the y-z plane, so the 
appropriate expression for the x component of the torque is 

dT x = r (x re R 2 sin 0 d0 d<{>) sin a sin <{> (43) 

From geometrical considerations, a = 0/2, and r = R{2( 1 - cos 0)} ^ , while x^ = (Sjiv^RJsin© 
[Bird, Stewart, and Lightfoot I960]. If the integration of Equation (1.5) is carried out over half a 
sphere (0 < <]> < n, 0 < 0 < 7t), then the x component of the torque is 
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and the corresponding force due to flow around half the sphere is SjiRjivj^. The force, F, on the 
hydrogen bond is twice the force due to flow around half the sphere, 

F = 37tR ^ v ~ (45) 

2 

which must equal the breaking force of the hydrogen bond if the molecule is denatured. 

From Hagler et al. , the weakest amide crystal bond (in the absence of the electrostatic 
contribution) has a breaking force of 8.3 x 10 13 N, which would require v^ = 9 cm/s to break the 
bond. This velocity is approximately two orders of magnitude greater than the free convection 
velocity one would estimate from the case of the vertical flat plate, and is certainly greater than the 


velocity attained in the systems of interest. Furthermore, this is a worst case scenario since the 
weakest possible bond was chosen and only one hydrogen bond was permitted. In reality, the 
electrostatic contribution would strengthen the bond and there are many internal hydrogen bonds. 

The analyses above indicate that viscous stresses due to natural convection are insufficient to 
disrupt crystal growth by stripping molecules from the surface, orienting molecules at the surface, 
or denaturing the protein as it approaches the surface. It is difficult to construct another mechanism 
by which the flow field can influence protein crystal growth except by altering the mass transfer 
rate to the crystal surface 

Effect of Mass Transfer Rate on Protein Crystal Growth 

The growth rate of tetragonal lysozyme crystals from solution (pH 4, salt concentration of 50 
mg/ml NaCl, 22° C) was measured by Pusey et al. [1986]. For small crystals (less than 70 |im in 
length), the growth rate was found to obey the relation 



where R is the crystal radius, Cj is the protein mass concentration at the interface, C s =1.7 mg/ml is 
the solubility mass concentration under these conditions, and k=1.46 x 10' 9 cm/s. No explanation 
of this unusual concentration dependence is given except that Schlichtkrull [1957] observed a 
similar behavior for insulin. The interfacial concentration calculated from their model of convective 
mass transfer was essentially the same as the bulk concentration. As we will show by means of a 
quasi-steady state analysis, the crystal sizes studied were too small for appreciable size effects to be 
evident and no reliable measurements have been made on larger crystals. Any discussion of 
growth rates of large crystals, then, is based on the extrapolation of small crystal measurements and 
should be considered speculative at best. 

The crystal growth rate given by Equation (2.1) must also equal the volume flux to the crystal 


surface: 
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where C x =725 mg/ml is the mass concentration of protein in the lysozyme crystal (corresponding 
to 50% solvent by volume [Matthews 1968, Bugg et al. 1984]). Note that the crystal shape has 
been approximated as spherical. The concentration gradient at the crystal surface can be expressed 
as 
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where Sh R is the Sherwood number based on the crystal radius and is given by the Ranz-Marshall 
correlation: 
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When the crystal grows by diffusion in a quasi-steady manner, the concentration gradient at the 
surface is given by 

C^-Ci 

= ” R— (50) 
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Equating the two expressions for crystal growth given in equations (46) and (47) and making use 
of equations (48) through (50) yields 


dC 

dr 


2 


[ 1 +0.3Sc 1/3 Gr 1/4 ] 


(51) 


i 


I 

+ 


Da 


C.-C. 
i s 


C -C. 


where Da=kR/D. Equation (51) can be solved for the interfacial concentration Cj as a function of 
size in order to determine the dependence of growth rate on crystal size. When g = 0 (no buoyancy 
driven flow). Equation (51) can be solved for Cj explicidy to yield 


C = 


•1 + 2 Da (Cy/Cg) + 1+ 4 Da (C x /C s )[(C oe - C s )/C $ ] 
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(52) 


The results of these calculations are given in Figure 6a, while the corresponding growth rate 
calculated according to Equation (46) is shown in Figure 6b. 

As shown in Figure 6a, the interfacial concentration is essentially equal to the bulk 
concentration over the entire size range studied by Pusey et al. whether or not natural convection is 
present The surface concentration and growth rate for crystals larger than 70 |im are stricdy 
extrapolations. If the crystal continues to grow at a nearly constant rate (broken lines), it would 
indicate that natural convection is sufficient to maintain the surface concentration at the bulk level so 
that crystal growth is entirely kinetically controlled. If convection is suppressed, a decrease in the 
growth rate from the small-size limit indicates that transport plays a role in controlling crystal 
growth as the crystal grows larger (solid lines). By the time a crystal growing from a 5%(w/v) 
solution of lysozyme has reached 1 mm in diameter, the growth rate has fallen to approximately 
25% of its initial value. Similar results can be seen for growth from l%(w/v) lysozyme solution. 

A rough indication of the relative importance of diffusion and kinetics can be obtained from the 
slope of the growth rate vs. crystal size curve (Figure 6b). The slope is zero when mass transfer is 
infinitely faster than interface kinetics, and the slope approaches -1 as diffusion begins to control 
the growth rate. From the curves in Figure 6b, it is clear that most crystal growth occurs in the 

transition region when both processes are important. 
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Quasi-steady models of this sort cannot predict the cessation of growth which has been 
reported [Kam et al. 1978, Feher and Kam 1985] because they do not include any mechanism 
which would reduce the ability of protein molecules to attach to the surface. The current data on the 
effect of mass transferrate on crystal growth are not sufficient to determine if mass transfer plays a 
role in limiting crystal size. If crystal defects occur because natural convection maintains the 
interfacial concentration at excessive levels, the simplest remedy may be to grow crystals from less 
concentrated solutions, in which case, the growth of high quality crystals may depend on the 
trade-off between growing large crystals and growing them quickly. 

Termination of Growth 

Growing sufficiently large crystals is often the limiting step in protein crystallography. 
Generally, protein crystals reach some terminal size beyond which they do not grow, even when 
they are transferred to a fresh solution of protein [C. Schutt 1986, personal conversation]. Several 
growth-ending mechanisms involving the effects of fluid flow were examined in Section 1 of this 
work, but they hardly exhaust the list of possible explanations. For example, the solution 
conditions may change during the course of crystal growth so that a transformation of the dissolved 
protein occurs. Cole et al. [ 1964] measured the pH of lysozyme solutions and found that during 
the course of crystallization, the pH increases by approximately 0 - 0.4 pH units from original 
values between pH 2.5 and pH 4.6. Ataka and Tanaka [1986], on the other hand, report a slight 
decrease in pH from solutions of pH 5 or higher which they attributed to absorption of carbon 
dioxide from air. Association of lysozyme molecules which might account for the lack of growth 
in the original solution has been observed in the range of pH 4.5 - 6.5 [Sophianopoulos and Van 
Holde 1964, Bruzzesi et al. 1965]. This mechanism, however, cannot account for the lack of 
growth when the crystals are placed in a fresh bath of protein unless the pH change also makes an 
irreversible change in the state of the protein on the crystal surface. 

The steady accumulation of errors suggested by Kam [Kam et al. 1978, Feher and Kam 1985] 


would explain why crystals stop growing, but the theory has not been verified. Experiments which 
study the relation between crystal size and defect concentration are necessary to confirm the validity 
of Kam’s hypothesis: if diffraction resolution does not increase with crystal size as expected, or if 
the resolution actually decreases, the accumulation of errors would be a suitable explanation. 
Finding a relation between defect concentration and crystal size, however, is not enough to explain 
how defects occur. To do this, a detailed study of the growing crystal is required so that local 
conditions such as interfacial concentrations of protein, precipitant, and hydrogen ion can be 
followed over the course of crystal growth. Such measurements, along with crystallographic 
studies of crystal ordering may provide some insight into the processes which terminate crystal 
growth. 

The mechanisms by which defects accumulate have not been studied because previous 
researchers worked only with small crystals (less then 70 pm in length) [Fiddis 1978, Pusey et al. 
1986]. Crystals this small would not produce strong size-dependent effects on the protein flux (see 
Ranz-Marshall correlation, Equation (2.4) ) and probably have a low concentration of surface 
defects. It follows that the kinetic expressions for the growth rate obtained from these 
measurements are incomplete because the role of defects has been neglected. Growth rate 
measurements taken over the full range of crystal sizes would provide information on the 
accumulation of defects and possibly on the evolution of any buoyancy-driven convection. The 
results of these experiments could then be used to evaluate the effect of natural convection on 
protein crystal growth and to explain why defects inhibit crystal growth. 
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FIGURE 1 - Quasi-equilibrium size distributions 

These distributions were calculated from the quasi-equilibrium approximation 
given by Equation (39) when the critical size is calculated from Equation (37) 
with K, = 0.065/%(w/v) and K« /K, = 35. 
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FIGURE 2 - Best-fit of trial size distribution 

Figure shows approximate best-fit of trial distribution 
obtained in the manner of Kam et al. and Kam and 
Feher. Solid curves were obtained from the quasi- 
equilibrium distribution for the values of the parameters; 
the linewidth was calculated from Equation (29). The 
open diamonds were calculated for the trial distribution. 
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FIGURE 3 - Free energy barrier to nucleation 

The maximum free energy required to form nuclei of 
critical size is calculated from Equation (38) for the 
reported values of = 0.065/%(w/v) and /K* =35. 
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FIGURE 4 - Effects of quasi-equilibrium approximation at high concentration 

(a) At concentrations where small changes in total protein concentration 
reduce j* by one, the amount of monomer in solution jumps 

discontinuously. • 

(b) The abrupt changes in monomer concentration affect the linewidth by 
making the system more monodisperse. 

All calculations were based on j* calculated from Equation (37), quasi- 
equilibrium concentrations calculated from Equation (39), and lhiewidths 
calculated from Equation (29). 
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hydrogen bond at A 
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FIGURE 5 - Denaturation of protein molecules by fluid shear. 

(a) Protein molecule in spherical (normal) conformation 
(b) Protein swinging open under influence of shear 
(c) Definition sketch showing coordinate system for calculations 
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FIGURE 6 - Surface concentration and growth rates of lysozyme crystals 

(a) Quasi-steady surface concentration of lysozyme when growth 
obeys the relation given by Equation (46). 

(b) Quasi-steady growth rate when growth obeys relation given . 

by Equation (46). r 

Growth rate constant reported by Pusey et al. to be k = 1.46x 10‘ 9 cm/s. 









